perm filename V2G.IN[TEX,DEK]1 blob sn#285512 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00021 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	folio 262 galley 1
C00016 00003	folio 265 galley 2
C00029 00004	folio 269 galley 3
C00041 00005	folio 273 galley 4
C00059 00006	folio 276 galley 5
C00073 00007	folio 282 galley 6
C00088 00008	folio 285 galley 7 WARNING: Much of this tape unreadable!
C00099 00009	folio 291 galley 8
C00112 00010	folio 293 galley 9
C00121 00011	folio 301 galley 10
C00137 00012	folio 300 galley 11
C00154 00013	folio 306 galley 12
C00174 00014	folio 313 galley 13
C00185 00015	folio 317 galley 14
C00206 00016	folio 322 galley 15-16 WARNING: End of 15 and beginning of 16 were unreadable.
C00222 00017	folio 328 galley 17
C00236 00018	folio 330 galley 18
C00247 00019	folio 333 galley 19
C00258 00020	folio 336 galley 20
C00269 00021	folio 338 galley 21
C00286 ENDMK
C⊗;
folio 262 galley 1
    0  {U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f. 
    2  262!ch. 4!g. 1a|'{A12}|π!|4|4|4|4The |¬m|¬i|¬x 
    7  computer assumes that its ⊗oating-point numbers 
   13  have the form|'{A9}|¬N!|9|εe!|9f!|9f!|9f!|9f!|9.|J!(4)|;
   17  {A9}|πHere we have base |εb |πexcess |εq |π⊗oating-point 
   25  notation with four ``digits'' of precision, where 
   32  |εb |πis the byte size (e.g., |εb|4α=↓|464 |πor 
   40  |εb|4α=↓|4100), |πand |εq |πis equal to |"l|f1|d32|)|εb|"L. 
   47  |πThe fraction part is |→|¬N|εffff, |πand |εe 
   54  |πis the exponent, which lies in the range |ε0|4|¬E|4e|4|¬W|
   62  4b. |πThis internal representation is typical 
   68  of the conventions in most existing computers, 
   75  although |εb |πis a much larger base than usual.|'
   84  !|4|4|4|4Most ⊗oating-point routines now in use 
   90  deal almost entirely with normalized numbers: 
   96  inputs to the routines are assumed to be normalized, 
  105  and the outputs are always normalized. Under 
  112  these conventions we lose the ability to represent 
  120  a few numbers of very small magnitude#for example, 
  128  the value (0,|4.00000001) can't be normalized 
  134  without producing a negative exponent#but we 
  140  gain in speed, uniformity, and the ability to 
  148  give relatively simple bounds on the relative 
  155  error in our computations. (Unnormalized ⊗oating-point 
  161  arithmetic is discussed in Section 4.2.2.)|'!|4|4|4|4Let 
  168  us now study the normalized ⊗oating-point arithmetic 
  175  in detail. At the same time we can consider the 
  185  construction of subroutines for these operations 
  191  , assuming that we have a computer without built-in 
  200  ⊗oating-point hardware.|'!|4|4|4|4Machine-language 
  203  subroutines for ⊗oating-point arithmetic are 
  208  usually written in a very machine-dependent manner, 
  215  using many of the wildest idiosyncrasies of the 
  223  computer at hand, and so we _nd that ⊗oating-point 
  232  addition subroutines for two di=erent machines 
  238  usually bear little super_cial resemblance to 
  244  each other. Yet a careful study of numerous subsroutines 
  253  for both binary and decimal computers reveals 
  260  that these programs actually have quite a lot 
  268  in common, and it is possible to discuss the 
  277  topics in a machine-independent way.|'!|4|4|4|4The 
  283  _rst (and by far the most di∃cult*3)) *?*?*?!|4|4|4|4The 
  291  _rst (and by far the most di∃cult*3) algorithm 
  299  we shall discuss in this section is a procedure 
  308  for ⊗oating-point addition:|'{A9}|ε(e|βu,|4f|βu)|4|↔V|4(e|βv
  311  ,|4f|βv)|4α=↓|4(e|βw,|4f|βw).|J!(6)|;{A9}Note|*/:|\ 
  313  Since ⊗oating-point arithmetic is inherently 
  318  approximate, not exact, we will use the symbols|'
  326  {A9}|↔V,!|↔B,!|↔N,!|↔n|;{A9}to denote ⊃oating-point 
  330  addition, subtraction, multiplication, and division, 
  335  respectively, in order to distinguish the approximate 
  342  operations from the true ones.|'{A12}|π!|4|4|4|4The 
  348  basic idea involved in ⊗oating-point addition 
  354  is fairly simple: Assuming that |εe|βu|4|¬R|4e|βv, 
  360  |πwe take |εe|βw|4α=↓|4e|βu, f|βw|4α=↓|4f|βu|4α+↓|4f|βv/b|ur
  363  e|βuα_↓e|βv|)|) |π(thereby aligning the radix 
  368  points for a meaningful addition), and normalize 
  375  the result. Several situations can arise which 
  382  make this process nontrivial, and the following 
  389  algorithm explains the method more precisely.|'
  395  |≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡A (|εFloating-point 
  398  addition).|9|4|πGiven base |εb, |πexcess |εq, 
  403  p-|πdigit, normalized ⊗oating-point numbers |εu|4α=↓|4(e|βu,
  407  |4f|βu) |πand |εv|4α=↓|4(e|βv,|4f|βv), |πthis 
  411  algorithm forms the sum |εw|4α=↓|4u|4|↔V|4v. 
  416  |πThe same procedure may be used for ⊗oating-point 
  424  subtraction, if |ε|→α_↓v |πis substituted for 
  430  |εv.|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡. |≡2|≡.|9|4Floating-point 
  433  addition.|;{A6}{H10L12M29}{I1.8H}|≡A|≡1|≡.|9[Unpack.] 
  435  Separate the exponent and fraction parts of the 
  443  representations of |εuu |π*?|≡A|≡1|≡.|9[Unpack.] 
  447  Separate the exponent and fraction parts of the 
  455  representations of |εu |πand |εv.|'{A3}|π|≡A|≡2|≡.|9[Assume 
  461  |εe|βu|4|¬R|4e|βv.] |πIf |εe|βu|4|¬W|4e|βv, |πinterchange 
  465  |εu |πand |εv. |π(In many cases, it is best to 
  475  combine step A2 with step A1 or with some of 
  485  the later steps.)|'{A3}|≡A|≡3|≡.|9[Set |εe|βw.] 
  490  |πSet |εe|βw|4|¬L|4e|βu.|'{A3}|π|≡A|≡4|≡.|9[Test 
  493  |εe|βu|4α_↓|4e|βv.] |πIf |εe|βu|4α_↓|4e|βv|4|¬R|4p|4α+↓|42 
  496  (|πlarge di=erence in exponents), set |εf|βw|4|¬L|4f|βu 
  502  |πand go to step A7. (Actually, since we are 
  511  assuming that |εu |πis normalized, we could terminate 
  519  the algorithm; but it is often useful to be able 
  529  to normalize a possibly unnormalized number by 
  536  adding zero to it.)|'{A3}|≡A|≡5|≡.|9[Scale right.] 
  542  Shift |εf|βv |πto the right |εe|βu|4α_↓|4e|βv 
  548  |πplaces; i.e., divide it by |εb|ure|βuα_↓e|βv|)|). 
  554  Note|*/: |\|πThis will be a shift of up to |εp|4α+↓|41 
  564  |πplaces, and the next step (which adds |εf|βu 
  572  |πto |εf|βv) |πthereby requires an accumulator 
  578  capable of holding |ε2p|4α+↓|41 |πbase |εb |πdigits 
  585  to the right of the radix point. If such a large 
  596  accumulator is not available, it is possible 
  603  to shorten the requirement to |εp|4α+↓|42 |πor 
  610  |εp|4α+↓|43 places if proper precautions are 
  616  taken; the details are given in exercise 5.|'
  624  {A3}|π|≡A|≡6|≡.|9[Add.] Set |εf|βw|4|¬L|4f|βu|4α+↓|4f|βv.|'
  627  {A3}|π|≡A|≡7|≡.|9[Normalize.] (At this point 
  631  (|εe|βw, f|βw) |πrepresents the sum of |εu |πand 
  639  |εv, |πbut |¬G|εf|βw|¬G |πmay have more than 
  646  |εp |πdigits, and it may be greater than utinty 
  655  or less than |ε1/b.) |πPerform Algorithm N below, 
  663  to normalize and round |ε(e|βw,|4f|βw) |πinto 
  669  the _nal anawer.|'{A12}{IC}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
  673  |≡N (|εNormalization).|9|4|πA ``raw exponent'' 
  677  |εe |πand a ``raw fraction'' |εf |πare converted 
  685  to normalized form, rounding if necessary to 
  692  |εp |πdigits. This algorithm assumes that |ε|¬G|1f|1|¬G|4|¬W
  698  |4b.|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡. |≡3|≡.|9|4Normalizati|π(`
  700  `fraction over⊗ow''), go to step N4. If |εf|4α=↓|40, 
  708  |πset |εe |πto its lowest possible value and 
  716  go to step N7.|'{A3}|≡N|≡5|≡.|9[Is |εf |πnormalized?] 
  723  If |ε|¬G|1f|1|¬G|4|¬R|41/b, |πgo to step N5.|'
  729  {A3}|≡N|≡3|≡.|9[Scale left.] Shift |εf |πto the 
  735  left by one digit position (i.e., multiply it 
  743  by |εb), |πand decrease |εe |πby 1. Return to 
  752  step N2.|'{A3}|≡N|≡4|≡.|9[Scale right.] Shift 
  757  |εf |πto the right by one digit position (i.e., 
  766  divide it by |εb), |πand increase |εe |πby 1.|'
  775  {A3}|≡N|≡5|≡.|9[Round.] Round |εf |πto |εp |πplaces. 
  781  (We take this to mean that |εf |πis changed to 
  791  the nearest multiple of |εb|gα_↓|gp. |πIt is 
  798  possible that (|εb|gpf)|πmod 1|4α=↓|4|f1|d32|) 
  802  so that there are |εtwo |πnearest multiples; 
  809  if |εb |πis even, we choose that one which makes 
  819  |εb|gpf|4α+↓|4|f1|d32|)b |πodd. Further discussion 
  823  of rounding appears in Section 4.2.2.) It is 
  831  important to note that this rounding operation 
  838  can make |ε|¬G|1f|1|¬G|4α=↓|41 (|π``rounding 
  842  over⊗ow''); in such a case, return to step N4.|'
  851  {A3}|≡N|≡6|≡.|9[Check |εe.] |πIf |εe |πis too 
  857  large, i.e., larger than its allowed range, an 
  865  |εexponent over⊃ow |πcondition is sensed. If 
  871  |εe |πis too small, an |εexponent under⊃ow |πcondition 
  879  is sensed. (See the discussion below; since the 
  887  result cannot be expressed as a normalized ⊗oating-point 
  895  number in the required range, special action 
  902  is necessary.)|'{A3}|≡N|≡7|≡.|9[Pack.] Put |εe 
  907  |πand |εf |πtogether into the desired output 
folio 265 galley 2
  914  {U0}{H10L12M29}91582#Computer Programming!(Knuth/A.-W.)!f. 
  916  265!ch. 4!g. 2a|'{A12}!|4|4|4|4The following 
  921  |¬m|¬i|¬x subroutines, for addition and subtraction 
  927  of numbers having the form (4), show how Algorithms 
  936  A and N can be expressed as computer programs. 
  945  The subroutines below are designed to take one 
  953  input |εu |πfrom symbolic location |¬a|¬c|¬c, 
  959  and the other input |εv |πcomes from register 
  967  A upon entrance to the subroutine. The output 
  975  |εw |πappears both in register A and location 
  983  |¬a|¬c|¬c. Thus, a _xed-point coding sequence|'
  989  {A9}{H9}|¬l|¬d|¬a|9|¬a;!|¬a|¬d|¬d|9|¬b;!|¬s|¬u|¬b|9|¬c;!|¬s|
  989  ¬t|¬a|9|¬d;|J!(7)|;{A9}{H10}would correspond 
  992  to the ⊗oating-point coding sequence|'{A9}{H9}|¬l|¬d|¬a|9|¬a
  997  |¬,!|¬s|¬t|¬a|9|¬a|¬c|¬c;!|¬l|¬d|¬a|9|¬b|¬,!|¬j|¬m|¬p|9|¬f|¬
  997  a|¬d|¬d;!|¬l|¬d|¬a|9|¬c,!|¬j|¬m|¬p|9|¬f|¬s|¬u|¬b;!|¬s|¬t|¬a|
  997  9|¬d.|J!(8)|;{A12}{H10}|≡P|≡r|≡o|≡g|≡r|≡a|≡m 
  999  |≡A (|εAddition, subtraction, and normalization).|9|4|πThe 
 1004  following program is a subroutine for Algorithm 
 1011  A, and it is also designed so that the normalization 
 1021  portion can be used by other subroutines which 
 1029  appear later in this section. In this program 
 1037  and in many other programs throughout this chapter, 
 1045  |¬o|¬f|¬l|¬o stands for a subroutine which prints 
 1052  out a message to the e=ect that |¬m|¬i|¬x's over⊗ow 
 1061  toggle was unexpectedly found to be ``on.'' The 
 1069  byle size |εb |πis assumed to be a multiple of 
 1079  4. The normalization routine |¬n|¬o|¬r|¬m assumes 
 1085  that rI2|4α=↓|4|εe |πand rAX|4α=↓|4|εf, |πwhere 
 1090  rA|4α=↓|40 implies rX|4α=↓|40 and rI2|4|¬W|4|εb.|'
 1095  {A12}{H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!!!|∂!!!!!!!!!!!!!!!!
 1095  !!!|4|4|4|∂|E|;|>|*/|↔c|↔O|\|'|π|¬e|¬x|¬p|'|¬e|¬q|¬u|'
 1100  |¬1|¬.|¬1|'De_nition of exponent _eld.|'>|>|ε|*/|↔c|↔P|\|'
 1108  |π|¬f|¬s|¬u|¬b|'|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'Floating-point 
 1112  subtraction subroutine:|'>|>|ε|*/|↔c|↔L|\|'|π|;
 1118  |¬l|¬d|¬a|¬n|'|¬t|¬e|¬m|¬p|'Change sign of operand.|'
 1124  |>|ε|*/|↔c|↔M|\|'|π|¬f|¬a|¬d|¬d|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'
 1129  Floating-point addition subroutine:|'>|>|ε|*/|↔c|↔C|\|'
 1135  |π|;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'Ensure over⊗ow is 
 1141  o=.|'>|>|ε|*/|↔c|↔o|\|'|π|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
 1148  |¬t|¬e|¬m|¬p|4|¬L|4|εv.|'>|>|*/|↔c|↔p|\|'|;|π|¬l|¬d|¬x|'
 1154  |¬a|¬c|¬c|'rX|4|¬L|4|εu.|'>|>|*/|↔c|↔l|\|'|;|π|¬c|¬m|¬p|¬a|'
 1161  |¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'|εSteps A|*/|↔O|\, A|*/|↔P|\, 
 1165  A|*/|↔L|\ are combined here|*/:|\|'>|>|*/|↔c|↔m|\|'
 1172  |;|π|¬j|¬g|¬e|'|¬1|¬f|'Jump if |εe|βv|4|¬R|4e|βu.|'
 1178  >|>|*/|↔O|↔c|\|'|;|π|¬s|¬t|¬x|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'
 1184  |¬f|¬u|4|¬L|4|¬N|εf|4f|4f|4f|40.|'>|>|*/|↔O|↔O|\|'
 1188  |;|π|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'rI2|4|¬L|4|εe|βw.|'
 1192  >|>|*/|↔O|↔P|\|'|;|π|¬s|¬t|¬a|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
 1198  |;>|>|ε|*/|↔O|↔L|\|'|;|π|¬l|¬d|¬1|¬n|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬
 1204  p|≤&|'rI1|4|¬L|4α_↓|εe|βv.|'>|>|*/|↔O|↔M|\|'|;
 1210  |π|¬j|¬m|¬p|'|¬4|¬f|'>|>|ε|*/|↔O|↔C|\|'|π|¬1|¬h|'
 1216  |¬s|¬t|¬a|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'|¬f|¬u|4|¬L|4|¬N|εf|4f|4f|
 1218  4f|40 (u, v |πinterchanged).|'>|>|ε|*/|↔O|↔o|\|'
 1225  |;|π|¬l|¬d|¬2|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|≤&|'rI2|4|¬L|4|εe|βw
 1228  .|'>|>|*/|↔O|↔p|\|'|;|π|¬s|¬t|¬x|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
 1235  >|>|ε|*/|↔O|↔l|\|'|;|¬l|¬d|¬i|¬n|'|π|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|
 1240  'rI1|4|¬L|4α_↓|εe|βv.|'>|>|*/|↔O|↔m|\|'|π|¬4|¬h|'
 1246  |¬i|¬n|¬c|¬1|'|¬0|¬,|4|¬2|'rI1|4|¬L|4|εe|βu|4α_↓|4e|βv. 
 1249  (|πStep A4 unnecessary.)|'>|>|ε|*/|↔P|↔c|\|'|π|¬5|¬h|'
 1256  |¬l|¬d|¬a|'|¬f|¬v|'|εA|*/|↔C|\. Scale right.|'
 1261  >|>|*/|↔P|↔O|'|\|;|¬e|¬n|¬t|¬x|'|¬0|'|πClear|4rX.|'
 1268  >|>|ε|*/|↔P|↔P|\|'|;|π|¬s|¬r|¬a|¬x|'|¬0|¬,|¬1|'
 1274  Shift right |εe|βu|4α_↓|4e|βv |πplaces.|'>|>|ε|*/|↔P|↔L|\|'
 1281  |π|¬6|¬h|'|¬a|¬d|¬d|'|¬f|¬u|'|εA|*/|↔o|\. Add.|'
 1286  >|>|*/|↔P|↔M|\|'|;|π|¬j|¬o|¬v|'|¬n|¬4|'|εA|*/|↔p|\. 
 1293  Normalize. |πJump if fraction over⊗ow.|'>|>|ε|*/|↔P|↔C|\|'
 1301  |;|¬j|¬x|¬z|'|¬n|¬o|¬r|¬m|'|πEasy case?|'>|>|ε|*/|↔P|↔o|\|'
 1309  |;|π|¬c|¬m|¬p|¬a|?|≤$|¬0|≤$|≤#1|¬.|¬1|≤&|'Is 
 1313  |εf |πnormalized?|'>|>|ε|*/|↔P|↔p|\|'|;|¬j|¬n|¬e|'
 1320  |¬n|¬5|'|πif so, round it.|'>|>|ε|*/|↔P|↔l|\|'
 1328  |;|π|¬s|¬r|¬c|'|¬5|'|¬GrX|¬G|4!|4|¬GrA|¬G.|'>
 1333  |>|ε|*/|↔P|↔m|\|'|;|π|¬d|¬e|¬c|¬x|'|¬1|'(rX is 
 1340  positive.)|'>|>|ε|*/|↔L|↔c|\|'|;|π|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
 1347  (Operands had opposite signs,|'>|>|ε|*/|↔L|↔O|\|'
 1354  |;|π|¬s|¬t|¬a|'|¬h|¬a|¬l|¬f|≤#|¬0|¬.|¬0|≤&|'registers 
 1358  must be adjusted|'>|>|ε|*/|↔L|↔P|\|'|;|¬l|¬d|¬a|¬n|'
 1366  |¬t|¬e|¬m|¬p|'|πbefore rounding and normalization.)|'
 1371  >|>|ε|*/|↔L|↔L|\|'|;|π|¬a|¬d|¬d|'|¬h|¬a|¬l|¬f|'
 1377  >|>|ε|*/|↔L|↔M|\|'|;|π|¬a|¬d|¬d|'|¬h|¬a|¬l|¬f|'
 1383  complement least signi_cant half.|'>|>|ε|*/|↔L|↔C|\|'
 1390  |;|π|¬s|¬r|¬c|'|¬4|'Jump into normalization routine.|'
 1397  >|>|ε|*/|↔L|↔o|\|'|;|π|¬j|¬m|¬p|'|¬n|¬3|¬a|'>|>
 1405  |ε|*/|↔L|↔p|\|'|π|¬h|¬a|¬l|¬f|'|¬c|¬o|¬n|'|¬1//|>
 1409  |ε|*/|↔L|↔p|\|'|π|¬h|¬a|¬l|¬f|'|¬c|¬o|¬n|'|¬1//|¬2|'
 1413  One half the word size. (Sign varies)|'>|>|ε|*/|↔L|↔l|\|'
 1423  |π|¬f|¬u|'|¬c|¬o|¬n|'|¬0|'Fraction part |εf|βu.|'
 1429  >|>|*/|↔L|↔m|\|'|π|¬f|¬v|'|¬c|¬o|¬n|'|¬0|'Fraction 
 1436  part |εf|βv.|'>|>|*/|↔M|↔c|\|'|π|¬n|¬o|¬r|¬m|'
 1442  |¬j|¬a|¬z|'|¬z|¬r|¬o|'|εN|*/|↔O|\. Test f.|'>|>
 1449  |*/|↔M|↔O|\|'|π|¬n|¬2|'|¬c|¬m|¬p|¬a|?|≤$|¬0|≤$|≤#|¬1|¬.|¬1|≤&
 1452  |'|εN|*/|↔P|\. Is f normalized|*/?|\|'>|>|*/|↔M|↔P|\|'
 1460  |;|π|¬j|¬n|¬e|'|¬n|≡5|'if leading byte nonzero, 
 1467  to N5.|'>|>|ε|*/|↔M|↔L|\|'|π|¬n|¬3|'|¬s|¬l|¬a|¬x|'
 1474  |¬1|'|εN|*/|↔L|\. Scale left.|'|>|*/|↔M|↔M|\|'|π|¬n|¬3|¬a|'
 1481  |¬d|¬e|¬c|¬2|'|¬1|'Decrease |εe |πby 1.|'>|>|ε|*/|↔M|↔C|\|'
 1490  |;|π|¬j|¬m|¬p|'|¬n|¬2|'Return to N2.|'>|>|ε|*/|↔M|↔o|\|'
 1499  |π|¬n|¬4|'|¬e|¬n|¬t|¬x|'|¬1|'|εN|*/|↔M|\. Scale 
 1504  right.|'>|>|*/|↔M|↔p|\|'|;|π|¬s|¬r|¬c|'|¬1|'shift 
 1512  right, insert ``1'' with proper sign.|'>|>|ε|*/|↔M|↔l|\|'
 1521  |;|π|¬i|¬n|¬c|¬2|'|¬1|'Increase |εe |πby 1.|'
 1528  >|>|ε|*/|↔M|↔m|\|'|π|¬n|¬5|'|¬c|¬m|¬p|¬a|?|≤$|¬b|¬y|¬t|¬e/|¬2
 1533  |≤$|≤#|¬5|¬.|¬5|≤&|'|εN|*/|↔C|\. Round.|'>|>|*/|↔C|↔c|\|'
 1539  |;|π|¬j|¬l|'|¬n|¬6|'is |¬Gtail|¬G|4|¬W|4|f1|d32|)|εb?|'
 1544  >|>|*/|↔C|↔O|\|'|;|π|¬j|¬g|'|¬5|¬f|'>|>|ε|*/|↔C|↔P|\|'
 1553  |;|π|¬j|¬x|¬n|¬z|'|¬5|¬f|'Is |¬Gtail|¬G|4|¬Q|4|f1|d32|)|εb?|
 1557  '>|>|*/|↔C|↔L|\|'|;|π|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
 1564  |¬Gtail|¬G|4α=↓|4|f1|d32|)|εb; round to odd.|'
 1568  >|>|ε|*/|↔C|↔M|\|'|;|π|¬l|¬d|¬x|'|¬t|¬e|¬m|¬p|≤#|¬4|¬.|¬4|≤&|
 1573  '>|>|ε|*/|↔C|↔C|'|\|;|π|¬j|¬x|¬o|'|¬n|¬6|'To N6 
 1582  if rX is odd.|'>|>|ε|*/|↔C|↔o|\|'|π|¬5|¬h|'|¬s|¬t|¬a|'
 1591  α/↓|≤%|¬1|≤#|¬0|¬.|¬0|≤&|'Store sign of rA.|'
 1596  >|>|ε|*/|↔C|↔p|\|'|;|π|¬i|¬n|¬c|¬a|'|¬b|¬y|¬t|¬e|'
 1602  Add |εb|gα_↓|g4 |πto |ε|¬G|1f|1|¬G. |π(Sign varies)|'
 1608  >|>|ε|*/|↔C|↔l|\|'|;|π|¬j|¬o|¬v|'|¬n|¬4|'Check 
 1615  for rounding over⊗ow.|'>|>|ε|*/|↔C|↔m|\|'|¬n|¬6|'
 1622  |¬j|¬2|¬n|'|¬e|¬x|¬p|¬u|¬n|'N|*/|↔o|\. Check e. 
 1627  |πUnder⊗ow if |εe|4|¬W|40.|'>|>|ε|*/|↔o|↔c|\|'
 1633  |¬n|¬7|'|¬e|¬n|¬t|¬x|'|¬0|¬,|¬2|'N|*/|↔p|\. Pack. 
 1638  |πrX|4|¬L|4|εe.|'>|>|ε|*/|↔o|↔O|\|'|;|π|¬s|¬r|¬c|'
 1644  |¬1|'>|>|ε|*/|↔o|↔P|\|'|π|¬z|¬r|¬o|'|¬d|¬e|¬c|¬2|'
 1650  |¬b|¬y|¬t|¬e|'rI2|4|¬L|4|εe|4α_↓|4b.|'>|>|ε|*/|↔o|↔L|\|'
 1655  |π|¬8|¬h|'|¬s|¬t|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔o|↔M|\|'
 1661  |π|¬e|¬x|¬i|¬t|¬f|'|¬j|¬2|¬n|'α/↓|'Exit, unless 
 1666  |εe|4|¬R|4b.|'>|>|ε|*/|↔o|↔C|\|'|π|¬e|¬x|¬p|¬o|¬v|'
 1671  |¬h|¬l|¬t|'|¬2|'Exponent over⊗ow detected|'>|>
 1678  |ε|*/|↔o|↔o|\|'|π|¬e|¬x|¬p|¬u|¬n|'|¬h|¬l|¬t|'|¬1|'
 1682  Exponent under⊗ow detected|'>|>|ε|*/|↔o|↔p|\|'
 1688  |¬a|¬c|¬c|'|¬c|¬o|¬n|'|¬0|'|πFloating-point accumulator|'
 1693  >|Httt*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.
folio 269 galley 3
 1695  -W.)!f. 269!ch. g. 3a|'{A12}!|4|4|4|4The rather 
 1701  long section of code from lines 25 to 37 is needed 
 1712  because |¬m|¬i|¬x has only a 5-byte accumulator 
 1719  for adding signed numbers while in general |ε2p|4α+↓|41|4α=↓
 1726  |49 |πplaces of accuracy are required by Algorithm 
 1734  A. The program could be shortened to about half 
 1743  its present length if we were willing to sacri_ce 
 1752  a little bit of its accuracy, but we shall see 
 1762  in the next section that full accuracy is important. 
 1771  Line 55 uses a nonstandard |¬m|¬i|¬x instruction 
 1778  de_ned in Section 4.5.2 The running time for 
 1786  ⊗oating-point addition and subtraction depends 
 1791  on several factors which are analyzed in Section 
 1799  4.2.4.|'!|44|4*?*?*?*?!|4|4|4|4Now let us consider 
 1804  multiplication and division, which are simpler 
 1810  than addition, and which are some hat similar 
 1818  to each other.|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 1822  |≡M (|εFloating-point multiplication or division).|9|4|πGive
 1826  n base |εb, |πexcess |εq, p-|πdigit, normalized 
 1833  ⊗oating-point numbers |εu|4α=↓|4(e|βu,|4f|βu) 
 1836  |πand |εv|4α=↓|4(e|βv,|4f|βv), |πthis algorithm 
 1840  forms the product |εw|4α=↓|4u|4|↔N|4v |πor the 
 1846  quotient |εw|4α=↓|4u|4|↔n|4v.|'{A3}{I1.9H}|π|≡M|≡1|≡.|9[Unpa
 1848  ck.] Separate the exponent and fraction parts 
 1855  of the representations of |εu |πand |εv. |π(Sometimes 
 1863  it is convenient, but not necessary, to test 
 1871  the operands for zero during this step.)|'{A3}|≡M|≡2|≡.|9[Op
 1878  erate.] Set|'|h|ε|∂e|βw|4|¬L|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α+↓|4
 1880  1,!!|∂f|βw|4|¬L|4(b|gα_↓|g1f|βu)/f|βv!!|∂|πfor|4|1multiplica
 1880  tion;|E|n|;{A8}(9)|E|?{B8}|L|εe|βw|4|¬L|4e|βu|4α+↓|4e|βv|4α_
 1882  ↓|4q,|Lf|βw|4|¬L|4f|βuf|βv|L|πfor|4|1multiplication;>
 1883  |E|'{A4}|ε|Le|βw|4|¬L|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α+↓|41,|Lf|β
 1884  w|4|¬L|4(b|gα_↓|g1f|βu)/f|βv|L|πfor|4|1division.>
 1885  |E|'{A12}!!|1|1(Since the input numbers are assumed 
 1892  to be normalized, it follows that either |εf|βw|4α=↓|40, 
 1900  |πor |ε1/b|g2|4|¬E|4|¬G|1f|βw|¬G|4|¬W|41, |πor 
 1903  a division-by-zero error has occurred.) If necessary, 
 1910  the representation of |εf|βw |πmay be reduced 
 1917  to |εpα+↓2 |πor |εpα+↓3 |πdigits at this point, 
 1925  as in exercise 5.|'|≡M|≡3|≡.|9[Normalize.]|9|4Perform 
 1930  Algorithm N on |ε(e|βw,|4f|βw) |πto normalize, 
 1936  round, and pack the result. (|εNote|*/: |\|πNormalization 
 1943  is simpler in this case, since scaling left occurs 
 1952  at most once, and rounding over⊗ow cannot occur 
 1960  after division.)|'{A12}{IC}!|4|4|4|4The following 
 1964  |¬m|¬i|¬x subroutines, which use the same conventions 
 1971  as Program A above, illustrate the machine considerations 
 1979  necessary in connection with Algorithm M.|'{A12}|≡P|≡r|≡o|≡g
 1985  |≡r|≡a|≡m |≡M (|εMultiplication and division).|'
 1990  {H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!|∂!!!!!!!!!!!!!!!!!!!!!|4
 1990  |4|4|∂|E|;|>|ε|*/|↔c|↔O|\|'|¬q|'|¬e|¬q|¬u|'|¬b|¬y|¬t|¬e/|¬2|'
 1996  |εQ |πis half the byte size.|'>|>|ε|*/|↔c|↔P|\|'
 2005  |¬f|¬m|¬u|¬l|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'|πFloating-point 
 2009  multiplication subroutine:|'>|>|ε|*/|↔c|↔L|\|'
 2014  |;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'|πEnsure over⊗ow is 
 2020  o=.|'>|>|ε|*/|↔c|↔M|\|'|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
 2027  |¬t|¬e|¬m|¬p|4|¬L|4v.|'>|>|ε|*/|↔c|↔C|\|'|;|¬l|¬d|¬x|'
 2033  |¬a|¬c|¬c|'|πrX|4|¬L|4|εu.|'>|>|ε|*/|↔c|↔o|\|'
 2038  |;|¬s|¬t|¬x|'|¬f|¬u|≤#|¬0|¬.|¬4|≤&|'|¬f|¬u|4|¬L|4|¬N|4f|4f|4
 2041  f|4f|40.|'>|>|ε|*/|↔c|↔p|\|'|;|¬l|¬d|¬1|'|¬t|¬e|¬m|¬p|≤#|¬e|¬
 2047  x|¬p|≤&|'>|>|*/|↔c|↔l|\|'|;|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤
 2053  &|'>|>|*/|↔c|↔m|\|'|;|¬i|¬n|¬c|¬2|'|→|≤*/|¬q|¬,|¬1|'
 2060  |πrI2|4|¬L|4|εe|βu|4α+↓|4e|βv|4α_↓|4q.|'>|>|*/|↔O|↔c|\|'
 2064  |;|¬s|¬l|¬a|'|¬1|'>|>|*/|↔O|↔O|\|'|;|¬m|¬u|¬l|'
 2072  |¬f|¬u|'|πMultiply |εf|βu |πtimes |εf|βv.|'>|>
 2079  |*/|↔O|↔P|\|'|;|¬j|¬m|¬p|'|¬n|¬o|¬r|¬m|'|πNormalize, 
 2084  round, and exit.|'>|>|ε|*/|↔O|↔L|\|'α/↓|'>|>|*/|↔O|↔M|\|'
 2094  |¬f|¬d|¬i|¬v|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'|πFloating-point 
 2098  division subroutine:|'>|>|ε|*/|↔O|↔C|\|'|;|¬j|¬o|¬v|'
 2105  |¬o|¬f|¬l|¬o|'|πEnsure over⊗ow is o=.|'>|>|ε|*/|↔O|↔o|\|'
 2113  |;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'|¬t|¬e|¬m|¬p|4|¬L|4|εv.|'
 2117  >|>|ε|*/|↔O|↔p|\|'|;|¬s|¬t|¬a|'|¬f|¬v|≤#|¬0|¬.|¬4|≤&|'
 2123  |¬f|¬v|4|¬L|4|¬N|4f|4f|4f|4f|40.|'>|>|ε|*/|↔O|↔l|\|'
 2127  |;|¬l|¬d|¬1|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|≤&|'>|>
 2132  |ε|*/|↔O|↔m|\|'|;|¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|≤&|'
 2136  >|>|ε|*/|↔P|↔c|\|'|;|¬d|¬e|¬c|¬2|'|→|≤*/|¬q|¬,|¬1|'
 2142  |πrI2|4|¬L|4|εe|βu|4α_↓|4e|βv|4α+↓|4q.|'>|>|ε|*/|↔P|↔O|\|'
 2146  |;|¬e|¬n|¬t|¬x|'|¬0|'>|>|ε|*/|↔P|↔P|\|'|;|¬l|¬d|¬a|'
 2154  |¬a|¬c|¬c|'>|>|ε|*/|↔P|↔L|\|'|;|¬s|¬l|¬a|'|¬1|'
 2161  |πrA|4|¬L|4|εf|βu.|'>|>|ε|*/|↔P|↔M|\|'|;|¬c|¬m|¬p|¬a|'
 2167  |¬f|¬v|≤#|¬1|¬.|¬5|≤&|'>|>|ε|*/|↔P|↔C|\|'|;|¬j|¬l|'
 2173  α/↓|→|≤%|¬3|'|πJump if |¬G|1f|βu|¬G|4|¬W|4|¬G|1f|βv|¬G.|'
 2177  >|>|ε|*/|↔P|↔o|\|'|;|¬s|¬r|¬a|'|¬1|'|πOtherwise, 
 2184  scle |εf|βu |πright|'>|>|ε|*/|↔P|↔p|\|'|;|¬i|¬n|¬c|¬2|'
 2192  |¬1|'|π!!and increase rI2 by 1.|'>|>|ε|*/|↔P|↔l|\|'
 2201  |;|¬d|¬i|¬v|'|¬f|¬v|'|πDivide|'>|>|ε|*/|↔P|↔m|\|'
 2208  |;|¬j|¬n|¬o|¬v|'|¬n|¬o|¬r|¬m|'|πNormalize, round, 
 2213  and exit.|'>|>|ε|*/|↔L|↔c|\|'|¬d|¬v|¬z|¬r|¬o|'
 2219  |¬h|¬l|¬t|'|¬3|'|πUnnormalized or zero divisor.|'
 2225  >{A12}{H10L12M29}!|4|4|4|4The most noteworthy 
 2229  feature of this program is the provision for 
 2237  division in lines 24<27, which is made in order 
 2246  to ensure enough accuracy to round the answer. 
 2254  If |ε|¬G|1f|βu|¬G|4|¬W|4|¬G|1f|βv|¬G, |πstraightforward 
 2257  application of Algorithm M would leave a result 
 2265  of the form ``|¬N0|4|εf|4f|4f|4f'' |πin register 
 2271  A, and this would not allow a proper rounding 
 2280  without a careful analysis of the remainder (which 
 2288  appears in register X). So the program computes 
 2296  |εf|βw|4|¬L|4f|βu/f|βv |πin this case, ensuring 
 2301  that |εf|βw |πis either zero or normalized in 
 2309  all cases; rounding can proceed with _ve signi_cant 
 2317  bytes, possibly testing whether the remainder 
 2323  is zero.|'!|4|4|4|4We occasionally need to convert 
 2330  between _xed- and ⊗oating-point representations. 
 2335  A ``_x-to-⊗oat'' routine is easily obtained with 
 2342  the help of the normalization algorithm above; 
 2349  for example, in |¬m|¬i|¬x, the following subroutine 
 2356  converts an integer to ⊗oating-point form:|'{A12}{H9L11M29}|
 2362  ∂!!!|∂!!!!|∂!!!!|∂!!!!!|∂!!!!!!!!!!!!!!!|∂|E|;
 2363  |>|ε|*/|↔c|↔O|\|'|¬f|¬l|¬o|¬t|'|¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬f|'
 2368  |πAssume that rA|4α=↓|4|εu, |πan integer.|'>|>
 2375  |ε|*/|↔c|↔P|'|\|;|¬j|¬o|¬v|'|¬o|¬f|¬l|¬o|'|πEnsure 
 2380  over⊗ow is o=.|'>|>|ε|*/|↔c|↔L|\|'|;|¬e|¬n|¬t|¬2|'
 2388  |¬q|≤%|¬5|'|πSet raw exponent.|'>|>|ε|*/|↔c|↔M|\|'
 2395  |;|¬n|¬t|¬x|'|¬0|'>|>|*/|↔c|↔C|\|'|;|¬j|¬m|¬p|'
 2403  |¬n|¬o|¬r|¬m|'|πN{U0}{H10L12M29}58320#Computer 
folio 273 galley 4
 2405  Programming!(Knuth/A.-W.)!f. 273!ch. 4!g. 4a|'
 2409  {A12}!|4|4|4|4The debugging of ⊗oating-point 
 2413  subroutines is usually a di∃cult job, since there 
 2421  are so many cases to consider. Here is a list 
 2431  of common pitfalls which often trap a programmer 
 2439  or machine designer who is preparing ⊗oating-point 
 2446  routines:|'{A3}!|4|4|4|41)|9|εLosing the sign.|4|9|πOn 
 2450  many machines (not |¬m|¬i|¬x), shift instructions 
 2456  between registers will a=ect the sign, and the 
 2464  shifting operations used in normalizing and scaling 
 2471  numbers must be carefully analyzed. The sign 
 2478  is also lost frequently when minus zero is present. 
 2487  (For example, Program A is careful to retain 
 2495  the sign of register A in lines 30<34. See also 
 2505  exercise 6.)|'!|4|4|4|42)|9|εFailure to treat 
 2510  exponent under⊗ow or over⊗ow properly.|4|9|πThe 
 2515  size of |εe|βw |πshould not be checked until 
 2523  |εafter |πthe rounding and normalization, because 
 2529  preliminary tests may give an erroneous indication. 
 2536  Exponent under⊗ow and over⊗ow can occur on ⊗oating-point 
 2544  addition and subtraction, not only during multiplication 
 2551  and division; and even though this is a rather 
 2560  rare occurrence, it must be tested each time. 
 2568  Enough information should be retained that meaningful 
 2575  corrective actions are possible after over⊗ow 
 2581  or under⊗ow has occurred.|'!|4|4|4|4It has unfortunately 
 2588  become customary in many instances to ignore 
 2595  exponent under⊗ow and simply to set under⊗owed 
 2602  results to zero with no indication of error. 
 2610  This causes a serious loss of accuracy in most 
 2619  cases (indeed, it is the loss of |εall |πthe 
 2628  signi_cant digits), and the assumptions underlying 
 2634  ⊗oating-point arithmetic have broken down, so 
 2640  the programmer really must be told when under⊗ow 
 2648  has occurred. Setting the result to zero is appropriate 
 2657  only in certain cases when the result is to be 
 2667  added to a signi_cantly larger quantity. When 
 2674  exponent under⊗ow is not detected, we _nd mysterious 
 2682  situations in which |ε(u|4|↔N|4v)|4|↔N|4w |πis 
 2687  zero, but |ε(u|4|↔N|4w)|4|↔N|4v |πis not, since 
 2693  |εu|4|↔N|4v |πresults in exponent under⊗ow but 
 2699  |ε(u|4|↔N|4w)|4|↔N|4v |πcan be calculated without 
 2704  any exponents falling out of range. Similarly, 
 2711  we can _nd positive numbers |εa, b, c, d, |πand 
 2721  |εy |πsuch that|'{A9}{A8}(11)|E|?{B8}|ε(a|4|↔N|4y|4|↔V|4b)|4
 2725  |↔n|4(c|4|↔N|4y|4|↔V|4d)|4|∂|¬V|4|f2|d33|),|;
 2726  {A4}| (a|4|↔V|4b|4|↔n|4y)|4|↔n|4(c|4|↔V|4d|4|↔n|4y)|4|Lα=↓|4
 2726  1>{A9}|πif exponent under⊗ow is not detected. 
 2733  (See exercise 9.) Even though ⊗oating-point routines 
 2740  are not precisely accurate, such a disparity 
 2747  as (11) is quite unexpected when |εa, b, c, d, 
 2757  |πand |εy |πare all |εpositive|π*3 Exponent under⊗ow 
 2764  is usually not anticipated by a programmer, so 
 2772  he needs to be told about it.|'!|4|4|4|43)|9|εInserted 
 2780  garbage.|4|9|πWhen scaling to the left it is 
 2787  important to keep from introducing anything but 
 2794  zeros at the right. For example, note the ``|¬e|¬n|¬t|¬x 
 2803  |¬0'' instruction in line 21 of Program A, and 
 2812  the all-too-easily-forgotten ``|¬e|¬n|¬t|¬x |¬0'' 
 2816  instruction in line 04 of the |¬f|¬l|¬o|¬t subroutine 
 2824  (10). (But it would be a mistake to clear register 
 2834  X after line 29 in the division subroutine*3|'
 2842  !|4|4|4|44)|9|εUnforeseen rounding over⊃ow.|4|9|πWhen 
 2845  a number like .999999997 is rounded to 8 digits, 
 2854  a carry will occur to the left of the decimal 
 2864  point, and the result must be scaled to the right. 
 2874  Many people have mistakenly concluded that rounding 
 2881  over⊗ow is impossible during multiplication, 
 2886  since they look at the maximum value of |ε|¬G|1f|βuf|βv|¬G 
 2895  |πwhich is |ε1|4α_↓|42b|gα_↓|gp|4α+↓|4b|gα_↓|g2|gp, 
 2898  |πand this cannot round up to 1. The fallacy 
 2907  in this reasoning is exhibited in exercise 11. 
 2915  Curiously, the phenomenon of rounding over⊗ow 
 2921  |εis |πimpossible during ⊗oatin-point division 
 2926  (see exercise 12).|'|ε!|4|4|4|4(Note|*/: |\|πThere 
 2931  is one school of though which says it is harmless 
 2941  to ``round'' .999999997 to .99999999 instead 
 2947  of to 1.0000000, since this does not increase 
 2955  the worst-case bounds on relative error. Furthermore 
 2962  the ⊗oating-point number 1.0000000 may be said 
 2969  to represent all real values in the interval 
 2977  [1.0000000|4α_↓|45|4α⊗↓|410|gα_↓|g8, 1.0000000|4α+↓|45|4α⊗↓|
 2978  410|gα_↓|g8], while .99999999 represents all 
 2983  values in the much smaller interval (.99999999|4α_↓|45|4α⊗↓|
 2989  410|gα_↓|g9, .99999999|4α+↓|45|4α⊗↓|410|gα_↓|g9). 
 2991  Even though the latter interval does not contain 
 2999  the original number .999999997, each number of 
 3006  the second interval is contained in the _rst, 
 3014  so subsequent calculations with the second interval 
 3021  are no less accurate than with the _rst. But 
 3030  this argument is incompatible with the mathematical 
 3037  philosopny of ⊗oating-point arithmetic which 
 3042  is expressed in Section 4.2.2.)|'!|4|4|4|45)|9|εRounding 
 3048  before normalizing.|4|9|πInaccuracies are caused 
 3052  by premature rounding in the wrong digit position. 
 3060  This error is obvious when rounding is being 
 3068  done to the left of the appropriate position; 
 3076  and it is also dangerous in the less obvious 
 3085  cases where rounding _rst is done too far to 
 3094  the right, then is followed by rounding in the 
 3103  true position. For this reason it is a mistake 
 3112  to round during the ``scaling-right'' operation 
 3118  in step A5, except as prescribed in exercise 
 3126  5. (Yet the special case of roundkcng*?prescribed 
 3133  in exercise 5. (Yet the special case of rounding 
 3142  in step N5, then rounding again after rounding 
 3150  over⊗ow has occurred, is harmless, because rounding 
 3157  over⊗ow has occurred, is harmless, because rounding 
 3164  over⊗ow always yields |→|¬N1.0000000 and this 
 3170  is una=ected by the subsequent rounding process.)|'
 3177  !|4|4|4|46)|9|εFailure to retain enough precision 
 3182  in intermediate calculations.|4|9|πDetailed analyses 
 3186  of the accuracy of ⊗oating-point arithmetic, 
 3192  made in the next section, suggest strongly that 
 3200  normalizing ⊗oating-point routines should always 
 3205  deliver a properly rounded result to the maximum 
 3213  possible accuracy. There should be no exceptions 
 3220  to this dictum, even in cases which occur with 
 3229  extremely low probability; the appropriate number 
 3235  of signi_cant digits should be retained throughout 
 3242  the computations, as stated in Algorithms A and 
 3250  M.|'{A12}|≡C|≡.|9|≡F|≡l|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡o|≡i|≡n|≡t 
 3252  |≡h|≡a|≡r|≡d|≡w|≡a|≡r|≡e|≡.|9|4Nearly every large 
 3255  computer intended for scienti_c calculations 
 3260  includes ⊗oating-point arithmetic as part of 
 3266  its repertoire of built-in operations. Unfortunately, 
 3272  the design of such hardware usually includes 
 3279  some anomalies which result in dismally poor 
 3286  behavior in certain circumstances, and we hope 
 3293  that future computer designers will pay more 
 3300  attention to providing the proper behavior than 
 3307  they have in the past. It cost only a little 
 3317  more to build the machine right, and considerations 
 3325  in the following section show that substantial 
 3332  bene_ts will be gained.|'!|4|4|4|4The |¬m|¬i|¬x 
 3338  computer, which is being used as an example of 
 3347  a ``typical'' machine in this series of books, 
 3355  has an optional ``⊗oating-point attachment'' 
 3360  (available at extra cost) which includes the 
 3367  following operations:|'|*2?|9|¬f|¬a|¬d|¬d|¬, |¬f|¬s|¬u|¬b|¬, 
 3371  |¬f|¬m|¬u|¬l|¬, |¬f|¬d|¬i|¬v|¬, |¬f|¬l|¬o|¬t|¬, 
 3374  |¬f|¬c|¬m|¬p (C|4α=↓|41, 2, 3, 4, 5, 56, respectively; 
 3382  F|4α=↓|46). The contents of rA after the operation 
 3390  ``|¬f|¬a|¬d|¬d |¬v'' are precisely the same as 
 3397  the contents of rA after the operations|'{A9}{H9L11}|h|∂|4|4
 3404  |4|4|4|4|4|4|1!|4|4|4|4|4|4|4|4|4|4|4|E|n|;|L|¬s|¬t|¬a!|¬a|¬
 3405  c|¬c>|L|¬l|¬d|¬a!|¬v>|L|¬j|¬m|¬p!|¬f|¬a|¬d|¬d>
 3408  {A9}{H10L12M29}where |¬f|¬a|¬d|¬d is the subroutine 
 3413  which appears earlier in this section, except 
 3420  that both operands are automatically normalized 
 3426  before entry to the subroutine, if they are not 
 3435  already in normalized form. (If exponent under⊗ow 
 3442  occurs during this pre-normalization, but not 
 3448  during the normalization of the answer, no under⊗ow 
 3456  is signalled.) Similar remarks apply to |¬f|¬s|¬u|¬b, 
 3463  |¬f|¬m|¬u|¬l, and |¬f|¬d|¬i|¬v. The contents 
 3468  of rA after the operation ``|¬f|¬l|¬o|¬t'' are 
 3475  the contents after ``|¬j|¬m|¬p |¬f|¬l|¬o|¬t'' 
 3480  in the subroutine (10) above.|'!|4|4|4|4The contents 
 3487  of rA are unchanged by the operation ``|¬f|¬c|¬m|¬p 
 3495  |¬v''; this instruction sets the comparison indicator 
 3502  to less, equal, or greater, depending on whether 
 3510  the contents of rA are ``de_nitely less than,'' 
 3518  ``approximately equal to,'' or ``de_nitely greater 
 3524  than'' |¬v; this subject is discussed in the 
 3532  next section, and the precise action is de_ned 
 3540  by the subroutine |¬f|¬c|¬m|¬p of exercise 4.2.2<17 
 3547  with |¬e|¬p|¬s|¬i|¬l|¬o|¬n in location 0.|'!|4|4|4|4No 
 3553  register other than rA is a=ected by any of the 
 3563  ⊗oating-point operations. If exponent over⊗ow 
 3568  or under⊗ow occurs, the over⊗ow toggle is turned 
 3576  on and the exponent of the answer is given modulo 
 3586  the byte size. Division by zero leaves unde_ned 
 3594  garbage in rA.|'!|4|4|4|4Sometimes it is helpful 
 3601  to use ⊗oating-point operators in a nonstandard 
 3608  manner. For example, if the operation |¬f|¬l|¬o|¬t 
 3615  had not been included as part of |¬m|¬i|¬x's 
 3623  ⊗oating-point attachment, we could easily achieve 
 3629  its e=ect on 4-byte numbers by writing|'{A9}{H9L11M29}{A28}(
 3636  12)|E|?{B28}|∂|4|4|4|4|4|4|4|4|4|4|4!|∂|4|4|4|4|4|4|4|4|4|4|
 3637  4!|∂|4|4|4|4|4|4|4|4|1|∂|E|;|>|¬f|¬l|¬o|¬t|'|¬s|¬t|¬j|'
 3641  |¬9|¬f|'>|>|;|¬s|¬l|¬a|'|¬1|'>|>|;|¬e|¬n|¬t|¬x|'
 3651  |¬q|≤%|¬4|'>|>|;|¬s|¬r|¬c|'|¬1|'>|>|;|¬f|¬a|¬d|¬d|'
 3661  |→|≤$|¬0|≤$|'>|>|¬9|¬h|'|¬j|¬m|¬p|'α/↓|'>{A12}{H10L12M29}Thi
 3668  s routine is not strictly equivalent to the |¬f|¬l|¬o|¬t 
 3677  operator, since it assumes that the 1|1|1:|1|11 
 3684  byte of rA is zero, and it destroys rX. The handling 
 3695  of more general situations is a little trickly, 
 3703  because rounding over⊗ow can occur even during 
 3710  a |¬f|¬l|¬o|¬t operation.|'|H*?{U0}{H10L12M29}58320#Computer 
folio 276 galley 5
 3714  Programming!(Knuth/A.-W.)!f. 276!ch. 4!g. 5a|'
 3718  {A12}!|4|4|4|4Similarly, if we want to round 
 3724  a number |εu |πfrom ⊗oating-point form to the 
 3732  nearest _xed-point integer, and if we know that 
 3740  number is nonnegative and will _t in at most 
 3749  three bytes, we may write|'{A9}{H9}|¬f|¬a|¬d|¬d!|¬f|¬u|¬d|¬g
 3754  |¬e|;{A9}{H10}where |¬f|¬u|¬d|¬g|¬e contains 
 3758  the constant|'{A12}{H9}|≤%!!|¬q|≤%|¬4!!|¬1!!|¬0!!|¬0!!|¬0!!;
 3760  |;{A9}{H10}the result in rA will be|'{A9}{H9}|≤%!!|¬q|≤%|¬4!
 3767  !|¬1!!round|4|ε(u)!!.|J!(13)|;{A9}{H10}|π!|4|4|4|4Some 
 3769  machines have ⊗oating-point hardware which suppresses 
 3775  automatic rounding and gives additional information 
 3781  about the less-signi_cant digits that would ordinarily 
 3788  be rounded o=. Such facilities are of use in 
 3797  extending the precision of ⊗oating-point calculations, 
 3803  but their implementation is usually subject to 
 3810  unfortunate anomalies, and unrounded results 
 3815  are generally undesirable in single-precision 
 3820  calculations.|'{A12}|≡D|≡.|9|≡H|≡i|≡s|≡t|≡o|≡r|≡y 
 3822  |≡a|≡n|≡d |≡b|≡i|≡b|≡l|≡i|≡o|≡g|≡r|≡a|≡p|≡h|≡y|≡.|4|9The 
 3824  origins of ⊗oating-point notation can be traced 
 3831  back to the Babylonians (c. 1800|4{H7}B{H10}.{H7}C{H10}.), 
 3837  who made use of radix 60 ⊗oating-point arithmetic 
 3845  but did not have a notation for the exponents. 
 3854  The appropriate exponent was always somehow ``understood'' 
 3861  by the man doing the calculations. At least one 
 3870  case has been found in which the wrong answer 
 3879  was given, because addition was performed with 
 3886  improper alignment of the operands, but such 
 3893  examples are very rare; see O. Neugebauer, |εThe 
 3901  Exact Sciences in Antiguity (|πPrinceton University 
 3907  Press, 1952), 26<27. Another early contribution 
 3913  to ⊗oating-point notation is due to the Greek 
 3921  mathematician Apollonius (3rd century {H7}B{H10}.{H7}C{H10}.
 3925  ), who apparently was the _rst to explain how 
 3934  to simplify multiplication by collecting powers 
 3940  of 10 separately from their coe∃cients, at least 
 3948  in simple cases. [For a discussion of Apollonius's 
 3956  method, see Pappus, |εMathematical Collections 
 3961  |π(4th century {H7}A{H10}.{H7}D{H10}.).] After 
 3965  the Babylonian civilization died out, the _rst 
 3972  signi_cant uses of ⊗oating-point notation for 
 3978  products and quotients did not emerge until much 
 3986  later, about the time logarithms were invented 
 3993  (1600) and shortly afterwards when Oughtred invented 
 4000  the slide rule (1630). The modern notation |ε``x|gn'' 
 4008  |πfor exponents was being introduced at about 
 4015  the same time; separate symbols for |εx |πsquared, 
 4023  |εx |πcubed, etc. had been in use before this.|'
 4032  !|4|4|4|4Floating-point arithmetic was incorporated 
 4036  into the design of some of the earliest computers. 
 4045  It was independently proposed by Leonardo Torres 
 4052  y Quevedo in Madrid, 1914; by Konrad Zuse in 
 4061  Berlin, 1936; and by George Stibitz in New Jersey, 
 4070  1939. Zuse's machines used a ⊗oating-binary representation 
 4077  which he called ``semi-legarithmic notation''; 
 4082  he also incorporated conventions for dealing 
 4088  with the special quantities ``|¬X'' and ``unde_ned.'' 
 4095  The _rst American computers to operate with ⊗oating-point 
 4103  arithmetic hardware were the Bell Laboratories' 
 4109  Model V and the Harvard Mark II, both of which 
 4119  were relay calculators designed in 1944. [See 
 4126  B. Randell, |εThe Origins of Digital Computers 
 4133  |π(Berlin: Springer, 1973), pp. 100, 155, 163<164, 
 4140  259<260; |εProc. Symp. on Large-Scale Digital 
 4146  Calculating Machinery |π(Harvard, 1947), 41<68, 
 4151  69<79; |εDatamation |≡1|≡3 |π(April 1967), 35<44, 
 4157  (May 1967), 45<49; |εZeitschrift f|=4ur angewandte 
 4163  Mathematik und Physik |≡1 (1950), 345<346.]|'
 4169  |π{H10L12M29}!|4|4|4|4The use of ⊗oating-binary 
 4173  arithmetic was seriously considered in 1944<1946 
 4179  by researchers at the Moore School in their plans 
 4188  for the _rst |εelectronic |πdigital computers, 
 4194  but it turned out to be much harder to implement 
 4204  ⊗oating-point circuitry with tubes than with 
 4210  relays. The group realized that scaling is a 
 4218  problem in programming, but felt that it is only 
 4227  a very small part of a total programming job 
 4236  which is usually worth the time and trouble it 
 4245  takes, since it tends to keep a programmer aware 
 4254  of the numerical accuracy he is getting. Furthermore, 
 4262  they argued that ⊗oating-point representation 
 4267  takes up valuable memory space, since the exponents 
 4275  must be stored, and that it becomes di∃cult to 
 4284  adapt ⊗oating-point arithmetic to multiple precision 
 4290  calculations. [See von Neumann's |εCollected 
 4295  Works |π|≡5 (New York: Macmillan, 1963), 43, 
 4302  73<74.] At this time, of course, they were designing 
 4311  the _rst stored-program computer and the second 
 4318  electronic computer, and their choice was to 
 4325  be either _xed point |εor |π⊗oating point, not 
 4333  both. They anticipated the coding of ⊗oating-binary 
 4340  routines, and in fact ``shift-left'' and ``shift-right'' 
 4347  instructions were put into their machine primarily 
 4354  to make such routines more e∃cient. The _rst 
 4362  machine to have both kinds of arithmetic in its 
 4371  hardware was apparently a computer developed 
 4377  at General Electric Company [see |εProc. |*/|↔P|\nd 
 4384  Symp. on Large-Scale Digital Calculating Machinery 
 4390  (|πHarvard, 1948), 65<69].|'!|4|4|4|4Floating-point 
 4394  subroutines and interpretive systems for early 
 4400  machines were coded by D. J. Wheeler and others, 
 4409  and the _rst publication of such routines was 
 4417  in |εThe Preparation of Programs for an Electronic 
 4425  Digital Computer |πby Wilkes, Wheeler, and Gill 
 4432  (Reading, Mass.: Addison-Wesley, 1951), subroutines 
 4437  A1<A11, pp. 35<37, 105<117. It is interesting 
 4444  to note that ⊗oating |εdecimal |πsubroutines 
 4450  are described here, although a binary computer 
 4457  was being used; in other words, the numbers were 
 4466  represented as 10|ε|gef, |πnot |ε2|gef, |πand 
 4472  therefore the scaling operations required multiplication 
 4478  or division by 10. On this particular machine 
 4486  such decimal scaling was about as easy as shifting, 
 4495  and the decimal approach greatly simpli_ed input-output 
 4502  conversions.|'!|4|4|4|4Most published references 
 4506  to the details of ⊗oating-point arithmetic routines 
 4513  are scattered in ``technical memorandums'' distributed 
 4519  by various computer manufacturers, but there 
 4525  have been occasional appearances of these routines 
 4532  in the open literature. Besides the reference 
 4539  above, see R. H. Stark and D. B. MacMillan, |εMath. 
 4549  Comp. |≡5 (1951), 86<92, |πwhere a plugboard-wired 
 4556  program is described; D. McCracken, |εDigital 
 4562  computer Programming |π(New York: Wiley, 1957), 
 4568  121<131; J. W. Carr III, |εCACM |≡2 |π(May 1959), 
 4577  10<15; W. G. Wadey, |εJACM |≡7 (1960), 129<139; 
 4585  |πD. E. Knuth, |εJACM |≡8 (1961), 119<128; |πO. 
 4593  Kesner, |εCACM |≡5 (1962), 269<271; |πF. P. Brooks 
 4601  and K. E. Iverson, |εAutomatic Data Processing 
 4608  |π(New York: Wiley, 1963), 184<199. |πFor a discussion 
 4616  of ⊗oating-point arithmetic from a computer designer's 
 4623  standpoint, see ``Floating-point operation'' 
 4627  by S. G. Campbell, in |εPlanning a computer System, 
 4636  |πed. by W. Buchholz (New York: McGraw-Hill, 
 4643  1962), 92<121. Additional references are given 
 4649  in Section 4.2.2.|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
 4653  {A12}{H9L11M29}|9|1|≡1|≡.|9|4[|ε|*/|↔O|↔c|\] |πHow 
 4655  would Avogadro's number and Planck's constant 
 4661  be represented in base 100, excess 50, four-digit 
 4669  ⊗oating-point notation? (This would be the representation 
 4676  used by |¬m|¬i|¬x, as in (4),{U0}{H10L12M29}58320#Computer 
folio 282 galley 6
 4682  Programming!(Knuth/A.-W.)!f. 282!ch. 4!g. 6a|'
 4686  {A12}{H9L11M29}|9|1|≡2|≡.|9|4[|ε|*/|↔O|↔P|\] |πAssume 
 4688  that the exponent |εe |πis constrained to lie 
 4696  in the range |ε0|4|¬E|4e|4|¬E|4E; |πwhat are 
 4702  the largest and smallest positive values which 
 4709  can be written as base |εb, |πexcess |εq, p-|πdigit 
 4718  ⊗oating-point numbers? What are the largest and 
 4725  smallest positive values which can be written 
 4732  as |εnormalized |π⊗oating-point numbers with 
 4737  these speci_cations?|'{A3}|9|1|≡3|≡.|9|4[|ε|*/|↔O|↔O|\] 
 4740  |πShow that if we are using normalized ⊗oating-binary 
 4748  arithmetic, there is a way to increase the precision 
 4757  slightly without loss of memory space: A |εp-|πbit 
 4765  fraction part can be represented using only |εp|4α_↓|41 
 4773  |πbit positions of a computer word, if the range 
 4782  of exponent values is decreased very slightly.|'
 4789  {A3}|9|1|≡4|≡.|9|4[|ε|*/|↔O|↔c|\] |πAssume that 
 4792  |εb|4α=↓|410, p|4α=↓|48. |πWhat result does Algorithm 
 4798  A give for (50, |→α+↓.98765432)|4|↔V|4(49, |→α+↓.33333333)? 
 4804  For (53, |→α_↓.99987654)|4|↔V|4(54, |→α+↓.10000000)? 
 4808  For (45, |→α_↓.50000001)|4|↔V|4(54, |→α+↓.10000000)?|'
 4812  {A3}|9|1|≡5|≡.|9|4[|ε|*/|↔P|↔M|\] |πLet us say 
 4816  that |εx|4|¬Z|4y |π(with respect to a given radix 
 4824  |εb) |πif |εx |πand |εy |πare real numbers satisfying 
 4833  the following conditions:|'{A9}|"l|εx/b|"L|4α=↓|4|"ly/b|"L;|
 4836  ;{A4}x|4|πmod|4|εb|4α=↓|40!!|πi=!!|εy|4|πmod|4|εb|4α=↓|40;|;
 4838  {A4}0|4|¬W|4x|4|πmod|4|εb|4|¬W|4|f1|d32|)b!!|πi=!!|ε0|4|¬W|4
 4838  y|4|πmod|4|εb|4|¬W|4|f1|d32|)b;|;{A4}x|4|πmod|4|εb|4α=↓|4|f1
 4839  |d32|)b!!|πi=!!|εy|4|πmod|4|εb|4α=↓|4|f1|d32|)b;|;
 4840  {A4}|f1|d32|)b|4|¬W|4x|4|πmod|4|εb|4|¬W|4b!!|πi=!!|ε|f1|d32|
 4840  )b|4|¬W|4y|4|πmod|4|εb|4|¬W|4b.|;{A9}|πProve 
 4842  that if |εf|βv |πis replaced by |εb|gα_↓|gp|gα_↓|g2F|βv 
 4849  |πbetween steps A5 and A6 of Algorithm A, where 
 4858  |εF|βv|4|¬Z|4b|gp|gα+↓|g2f|βv, |πthe result of 
 4862  that algorithm will be unchanged. (If |εF|βv 
 4869  |πis an integer and |εb |πis even, this operation 
 4878  essentially truncates |εf|βv |πto |εp|4α+↓|42 
 4883  |πplaces while remembering whether any nonzero 
 4889  digits have been dropped, thereby minimizing 
 4895  the length of register which is needed for the 
 4904  addition in step A6.)|'{A3}|9|1|≡6|≡.|9|4[|ε|*/|↔P|↔c|\|π] 
 4909  If the result of a |¬f|¬a|¬d|¬d instruction is 
 4917  zero, what will be the sign of rA, according 
 4926  to the de_nitions in this section?|'|9|1|≡7|≡.|9|4[|ε|*/|↔P|↔
 4932  p|\|π] Discuss ⊗oating-point arithmetic using 
 4937  balanced ternary notation.|'{A3}|9|1|≡8|≡.|9|4[|ε|*/|↔P|↔c|\|
 4940  π] Give examples of normalized eight-digit ⊗oating-decimal 
 4947  numbers |εu |πand |εv |πfor which addition yields 
 4955  (a) exponent under⊗ow, (b) exponent over⊗ow, 
 4961  assuming that exponents must satisfy |ε0|4|¬E|4e|4|¬W|4100.|
 4966  '{A3}|π|1|9|≡9|≡.|9|4[|εM|*/|↔P|↔M|\] |π(W. Kahan.) 
 4970  Assume that the occurrence of exponent under⊗ow 
 4977  causes the result to be replaced by zero, with 
 4986  no error indication given. Using excess zero, 
 4993  eight-digit ⊗oating decimal numbers with |εe 
 4999  |πin the range |→α_↓50|4||¬E|4|εe|4|¬W|450, |π_nd 
 5004  positive values of |εa, b, c, d, |πand |εy |πsuch 
 5014  that (11) holds.|'{A3}|≡1|≡0|≡.|9|4[|ε|*/|↔O|↔P|\] 
 5018  |πGive an example of normalized eight-digit ⊗oating-decimal 
 5025  numbers |εu |πand |εv |πfor which rounding over⊗ow 
 5033  occurs in addition.|'{A3}|≡1|≡1|≡.|9|4[|εM|*/|↔P|↔c|\] 
 5037  |πGive an example of normalized excess 50, eight-digit 
 5045  ⊗oating-decimal numbers |εu |πand |εv |πfor which 
 5052  rounding over⊗ow occurs in multiplication.|'{A3}|≡1|≡2|≡.|9|
 5057  4[|εM|*/|↔P|↔C|\] |πProve that rounding over⊗ow 
 5062  cannot occur during the normalization phase of 
 5069  ⊗oating-point division.|'{A3}|≡1|≡3|≡.|9|4[|ε|*/|↔L|↔c|\] 
 5072  |πWhen doing ``interval arithmetic'' we don't 
 5078  want to round the results of a ⊗oating-point 
 5086  computation; we want rather to implement operations 
 5093  such as ! and ! which give the tightest possible 
 5103  representable bounds on the true sum: |εu|4!|4v|4|¬E|4u|4α+↓
 5109  |4v|4|¬E|4u|4!|4v. |πHow should the algorithms 
 5114  of this section be modi_ed for such a purpose?|'
 5123  {A3}|≡1|≡8|≡.|9|4[|ε|*/|↔P|↔C|\] |πConsider a 
 5126  binary computer with 36-bit words, on which positive 
 5134  ⊗oating-binary numbers are represented as (0|εe|β1e|β2|4.|4.
 5139  |4.|4e|β8f|β1f|β2|4.|4.|4.|4f|β2|β7)|β2; |πhere 
 5141  (|εe|β1e|β2|4.|4.|4.|4e|β8)|β2 |πis an excess 
 5145  (10000000)|β2 exponent and (|εf|β1f|β2|4.|4.|4.|4f|β2|β7)|β2
 5148   |πis a 27-bit fraction. Negative ⊗oating-point 
 5155  numbers are represented by the |εtwo's complement 
 5162  |πof the corresponding positive representation 
 5167  (see Section 4.1). This, 1.5 is (|ε|*/|↔P|↔c|↔O|↔o|↔c|↔c|↔c|↔
 5173  c|↔c|↔c|↔c|↔c|\)|β8 |πwhile |→α_↓1.5 is (|ε|*/|↔C|↔p|↔o|↔P|↔c
 5177  |↔c|↔c|↔c|↔c|↔c|↔c|↔c|\)|β8; |πthe octal representations 
 5181  of 1.0 and |→α_↓1.0 |ε|*/|↔P|↔c|↔O|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c
 5185   |\|πand |ε|*/|↔C|↔p|↔o|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔|\, 
 5188  |πrespectively. Note that bit |εf|β1 |πof a normalized 
 5196  positive number is always 1, while it is almost 
 5205  always zero for negative numbers; the exceptional 
 5212  cases are representations of |→α_↓2|ε|gk.|'|π!!|1|1Suppose 
 5218  that the exact result of a ⊗oating-point operation 
 5226  has the octal code |ε|*/|↔C|↔p|↔P|¬G|↔p|↔M|↔c|↔c|↔c|↔c|↔c|↔c|
 5230  ↔c|\|¬G|*/|↔c|↔O|\; |πthis (negative) 33-bit fraction 
 5235  must be normalized and rounded to 27 bits. If 
 5244  we shift left until the leading fraction bit 
 5252  is zero, we get |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c
 5256  |↔c|\|¬G|*/|↔P|↔c|\, |πbut this rounds to the 
 5262  illegal value |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔
 5264  c|\; |πwe have over-normalized, since the correct 
 5271  answer is |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\.
 5273   |πOn the other hand if we star with the value 
 5284  |ε|*/|↔C|↔p|↔P|\|¬G|*/|↔p|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\|¬G|*/|↔c|↔C
 5284  |\|π and stop before over-normalizing it, we 
 5291  obtain |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|\|¬G|
 5292  */|↔C|↔c|\ |πwhich rounds to the unnormalized 
 5298  number |ε|*/|↔C|↔p|↔C|\|¬G|*/|↔M|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔O|\|π; 
 5300  subsequent normalization yields |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔
 5303  c|↔c|↔c|↔c|↔c|↔c|↔P|\ |πwhile the correct answer 
 5308  is |ε|*/|↔C|↔p|↔o|\|¬G|*/|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔c|↔O.|'
 5310  !|1|1Give a simple, correct rounding rule which 
 5317  resolves this dilemma on such a machine (without 
 5325  abandoning two's complement notation).|'{A3}|≡1|≡9|≡.|9|4[|ε
 5329  |*/|↔P|↔M|\|π] What is the running time for the 
 5337  |¬f|¬a|¬d|¬d subroutine in Program A, in terms 
 5344  of relevant characteristics of the data? What 
 5351  is the maximum running time, over all inputs 
 5359  that do not cause over⊗ow or under⊗ow?|'{A3}|≡1|≡4|≡.|9|4[|ε
 5366  |*/|↔P|↔C|\|π] Write a |¬m|¬i|¬x subroutine which 
 5372  begins with an arbitrary ⊗oating-point number 
 5378  in register A, not necessarily normalized, and 
 5385  which converts it to the nearest _xed-point integer 
 5393  (or determines that it is too large in absolute 
 5402  value to make such a conversion possible).|'{A3}|≡1|≡5|≡.|9|
 5409  4[|ε|*/|↔P|↔l|\|π] Write a |¬m|¬i|¬x subroutine, 
 5414  to be used in connection with the other subroutines 
 5423  of this section, that calculates |εu|9|πmod|91, 
 5429  that is, |εu|4α_↓|4|"lu|"L rounded to nearest 
 5435  ⊗oating-point number, given a ⊗oating-point number 
 5441  |εu. |πNote that when |εu |πis a very smalll 
 5450  negative number, |εu|9|πmod|91 will be rounded 
 5456  so that the result is unity (even though |εu 
 5465  |πmod 1 has been de_nedf tto that the result 
 5474  is unity (even though |εu |πmod 1 has been de_ned 
 5484  to be always |εless |πthan unity, as a real number).|'
 5494  {A3}|≡1|≡6|≡.|9|4[|ε|*/HM|↔P|↔O|\|π] (Robert L. 
 5497  Smith.) Design an algorithm to compute the real 
 5505  and imaginary parts of the complex number |ε(a|4α+↓|4bi)/(c|
 5512  4α+↓|4di), |πgiven real ⊗oating-point values 
 5517  |εa, b, c, |πand |εd. |πAvoid the computation 
 5525  of |εc|g2|4α+↓|4d|g2, |πsince it would cause 
 5531  ⊗oating-point over⊗ow even when |ε|¬Gc|¬C |πor 
 5537  |ε|¬Gd|¬G |πis approximately the square root 
 5543  of the maximum allowable ⊗oating-point value.|'
 5549  {A3}|≡1|≡7|≡.|9|4[|ε|*/|↔M|↔c|\|π] (John Cocke.) 
 5552  Explore the idea of extending the range of ⊗oating-point 
 5561  numbers by de_ning a single-word representation 
 5567  in which the precision of the fraction decreases 
 5575  as the magnitude of the exponent increases.|'
folio 285 galley 7 WARNING: Much of this tape unreadable!
 5582  |{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f. 
 5584  285!ch. 4!g. 7a|'{A18}|∨4|∨.|∨2|∨.|∨2|∨.|9|∨A|∨c|∨c|∨u|∨r|∨a
 5587  |∨c|∨y|9|∨o|∨f|9|∨F|∨l|∨o|∨a|∨t|∨i|∨n|∨g|∨-|∨P|∨o|∨i|∨n|∨t|9
 5587  |∨A|∨r|∨i|∨t|∨h|∨m|∨e|∨t|∨i|∨c|'{A6}{H10L12M29}Floating-poin
 5588  t computation is by nature inexact, and it is 
 5597  not di∃cult to misuse it so that the computed 
 5606  answers consist almost entirely of ``noise.'' 
 5612  One of the principal problems of numerical analysis 
 5620  is to determine how accurate the results of certain 
 5629  numerical methods will be; a ``credibility-gap'' 
 5635  problem is involved here: we don't know how much 
 5644  of the computer's answers to believe. Novice 
 5651  computer users solve this problem by implicitly 
 5658  trusting in the computer as an infallible authority; 
 5666  they tend to believe that all digits of a printed 
 5676  answer are signi_cant. Disillusioned computer 
 5681  users have just the opposite approach, they are 
 5689  constantly afraid their answers are almost meaningless. 
 5696  Many a serious mathematician has attempted to 
 5703  give rigourous analyses of a sequence of ⊗oating-point 
 5711  operations, but has fouhd the task to be so formidable 
 5721  he has tried to content himself with plausibility 
 5729  arguments instead.|'!|4|4|4|4A thorough examination 
 5734  of error analysis techniques is, of course, beyond 
 5742  the scope of this book, but we will in this section 
 5753  study some of the characteristics of ⊗oating-point 
 5760  arithmetic errors. Our goal is to discover how 
 5768  to perform ⊗oating-point arithmetic in such a 
 5775  way that reasonable analyses of error propagation 
 5782  are facilitated as much as possible.|'!|4|4|4|4A 
 5789  rough (but reasonably useful) way to express 
 5796  the behavior of ⊗oating-point arithmetic can 
 5802  be based on the concept of ``signi_cant _gures'' 
 5810  or |εrelative error. |πIf we are representing 
 5817  an exact real number |εx |πinside a computer 
 5825  by using the approximation |ε|=7x|4α=↓|4x(1|4α+↓|4|≤e), 
 5830  |πthe quantity |ε|≤e|4α=↓|4(|=7x|4α_↓|4x)/x |πis 
 5834  called the relative error of approximation. Roughly 
 5841  speaking, the operations of ⊗oating-point multiplication 
 5847  and division do not magnify the relative error 
 5855  by very much; but ⊗oating-point subtraction of 
 5862  nearly equal quantities (and ⊗oating-point addition, 
 5868  |εu|4|↔V|4v, |πwhere |εu |πis nearly equal to 
 5875  |→α_↓|εv) |πcan bery greatly increase the relative 
 5882  error. So we have a general rule of thumb, that 
 5892  a substantial loss of accuracy is expected from 
 5900  such additions and subtractions, but not from 
 5907  multiplications and divisions. On the other hand, 
 5914  the situation is somewhat paradoxical and needs 
 5921  to be understood properly, since ``bad'' additions 
 5928  and subtractions are performed with perfect accuracy*3 
 5935  (See exercise 25.)|'!|4|4|4|4One of the consequences 
 5942  of the possible unreliability of ⊗oating-point 
 5948  addition is that the associative law breaks down:|'
 5956  {A9}|ε(u|4|↔V|4v)|4|↔V|4w|4|=|↔6α=↓|4u|4|↔V|4(v|4|↔V|4w),!!|
 5956  πfor|4|1certain|4|1|εu,|4v,|4w.|J!(1)|;{A9}|πFor 
 5958  example,|'{A9}(11111113.|4|↔V|4|→α_↓11111111.)|4|↔V|47.51111
 5959  11|4α=↓|42.0000000|4|↔V|47.5111111>{A3}α=↓|49.5111111;>
 5961  {A4}11111113.|4|↔V|4(|→α_↓11111111.|4|↔V|47.5111111)|4α=↓|41
 5961  1111113.|4|↔V|4|→α_↓11111103.>{A3}α=↓|410.000000.>
 5963  {A9}{H10L12M29}*?*?*?ct operations α+↓, α_↓, α⊗↓, 
 5968  /.)|'!|4|4|4|4In view of the failure of the associative 
 5977  law, the comment of Mrs La Touchehwhich appears 
 5985  at the beginning of this chapter [taken from 
 5993  |εMath. Gazette |≡1|≡2 (1924), 95] |πmakes a 
 6000  good deal of sense with respect to ⊗oating-point 
 6008  arithmetic. Mathematical notations like ``|εa|β1|4α+↓|4a|β2|
 6012  4α+↓|4a|β3'' |πor ``|¬K|ur|)1|¬Ek|¬E|εn|)|4a|βk'' 
 6015  |πare inherently based upon the assumption of 
 6022  associativity, so a programmer must be especially 
 6029  careful that he does not implicitly assume the 
 6037  validity of the associative law.|'{A12}|≡A|≡.|9|≡A|≡n 
 6043  |≡a|≡x|≡i|≡o|≡m|≡a|≡t|≡i|≡c |≡a|≡p|≡p|≡r|≡o|≡a|≡c|≡h|≡.|9|4A
 6044  lthough the associative law is not valid, the 
 6052  commutative law|'{A9}|εu|4|↔V|4v|4α=↓|4v|4|↔V|4u|J!(2)|;
 6055  {A9}|πdoes hold, and this law can be a valuable 
 6064  conceptual asset in programming and in the analysis 
 6072  of programs. This example suggests that we would 
 6080  look for important laws which |εare |πsati_ed 
 6087  by |↔V, |↔B, |↔N, and |↔n; it is not unreasonable 
 6097  to say that |ε⊃oating-point routines should be 
 6104  designed to preserve as many of the ordinary 
 6112  mathematical laws as possible.|'|π!|4|4|4|4Let 
 6117  us now consider some of the other basic laws 
 6126  which are valid for normalized ⊗oating-point 
 6132  operations as described in the previous section. 
 6139  First we have|'{A9}|h|ε|→α_↓(u|4|↔V|4v)|4|∂|4α_↓u|4|↔V|4|→α_
 6142  ↓v;|E|n|;| u|4|↔B|4v|4|Lα=↓|4u|4|↔V|4|→α_↓v;|J!(3)>
 6144  {A4}|→α_↓(u|4|↔V|4v)|4α=↓|4α_↓u|4|↔V|4|→α_↓v;|J!(4)|;
 6145  {A4}u|4|↔V|4u|4|↔N|4v|4α=↓|4|πround(|εu|4α⊗↓|4v),!!u|4|↔n|4v
 6145  |4α=↓|4|πround(|εu/v),|;{A9}|πwhere round(|εx) 
 6148  |πdenotes the best ⊗oating-point approximation 
 6153  to |εx. |πWe have|'{A9}round(|→α_↓|εx)|4α=↓|4α_↓|πround|ε(x)
 6157  ,|J!(10)|;{A4}x|4|¬E|4y!!|πimplies!!round(|εx)|4|¬E|4|πround
 6158  (|εy),|J!(11)|;{A9}|πand these fundamental relations 
 6163  lead to imediate proofs of properties (2) through 
 6171  (8). We can also write down several more identities:|'
 6180  {A9}|εu|4|↔N|4v|4α=↓|4v|4|↔N|4u,!!(|→α_↓u)|4|↔N|4v|4α=↓|4α_↓
 6180  (u|4|↔N|4v),!!1|4|↔N|4v|4α=↓|4v;|;{A5}u|4|↔N|4v|4α=↓|40!!|πi
 6181  f|4and|4only|4if!!|εu|4α=↓|40!!|πor!!|εv|4α=↓|40;|;
 6182  {A4}(|→α_↓u)|4|↔n|4v|4α=↓|4u|4|↔n|4(|→α_↓v)|4α=↓|4α_↓(u|4|↔n
 6182  |4v),!!0|4|↔n|4v|4α=↓|40,|;{A4}u|4|↔n|41|4α=↓|4u,!!u|4|↔n|4u
 6183  |4α=↓|41;|;{A9}|πif |εu|4|¬E|4v |πand |εw|4|¬Q|40 
 6188  |πthen |εu|4|↔N|4w|4|¬E|4v|4|↔N|4w, u|4|↔n|4w|4|¬E|4v|4|↔n|4
 6190  w, |πand |εw|4|↔n|4u|4|¬R|4w|4|↔n|4v. |πIf |εu|4|↔V|4v|4α=↓|
 6194  4u|4α+↓|4v |πthen |ε(u|4|↔V|4v)|4|↔B|4v|4α=↓|4u, 
 6197  |πand if |εu|4|↔N|4v|4α=↓|4u|4α⊗↓|4|=|↔6α=↓|40 
 6200  |πthen |ε(u|4|↔N|4v)|4|↔n|4v|4α=↓|4u, |πetc. 
 6203  We see that a good deal of regularity is present 
 6213  in spite of the inexactness of the ⊗oating-point 
folio 291 galley 8
 6221  o *?*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.
 6223  )!f. 291!ch. 4!g. 8a|'{A12}!|4|4|4|4Several familiar 
 6229  rules of algebra are still, of course, conspicuously 
 6237  absent from the collection of identities above; 
 6244  the associative law for ⊗oating-point multiplication 
 6250  is not strictly true, as shown in exercise 3, 
 6259  and the distributive law between |↔N and |↔V 
 6267  can fail rather badly: Let |εu|4α=↓|420000.000, 
 6273  v|4α=↓|4|→α_↓6.0000000, |πand |εw|4α=↓|46.0000003; 
 6276  |πthen|'|h|ε(u|4|↔N|4v)|4|↔N|4(u|4|↔N|4w)|4|∂α=↓|499999.000|
 6277  4|↔N|4.0000000000000|4α=↓|4.0000000000|E|n|;{A9}| (u|4|↔N|4v
 6278  )|4|↔V|4(u|4|↔N|4w)|4|Lα=↓|4|→α_↓120000.00|4|↔V|4120000.01|4
 6278  α=↓|4.010000000>{A4}| u|4|↔N|4(v|4|↔V|4w)|4|Lα=↓|420000.000|
 6279  4|↔N|4.00000030000000|4α=↓|4.0060000000>{A9}|πso|'
 6281  {A9}u|4|↔N|4(v|4|↔V|4w)|4|=|↔6α=↓|4(u|4|↔N|4v)|4|↔V|4(u|4|↔N
 6281  |4w).|J!(12)|;{A9}|πOn the other hand we do have 
 6289  |εb|4|↔N|4(v|4|↔V|4w)|4α=↓|4(b|4|↔N|4v)|4|↔V|4(b|4|↔N|4w), 
 6290  |πsince|'{A9}round(|εbx)|4α=↓|4b|4|πround(|εx).|J!(13)|;
 6292  {A9}|π(Strictly speaking, the identities and 
 6297  inequalities we are considering in this section 
 6304  implicitly assume that no exponent under⊗ow or 
 6311  over⊗ow occurs. The function round(|εx) |πis 
 6317  unde_ned when |ε|¬Gx|¬G |πis too small or too 
 6325  large, and equations such as (13) hold only when 
 6334  both sides are de_ned.)|'!|4|4|4|4The failure 
 6340  of Cauchy's fundamental inequality |ε(x|=|β1|g2|4α+↓|4|¬O|4|
 6344  ¬O|4|¬O|4α+↓|4x|=|βn|g2)(y|=|β1|g2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|
 6344  4y|=|βn|g2)|4|¬R|4(x,y|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4x|βny|βn
 6344  )|g2 |πis another important example where traditional 
 6351  algebra breaks down in the presence of ⊗oating-point 
 6359  arithmetic; exercise 7 exhibits ⊗oating-binary 
 6364  numbers |εu |πand |εv |πsuch that |ε2(u|=|g|*2`|g2|4|↔V|4v|=|
 6370  g|*2`|g2)|4|¬W|4(u|4|↔V|4v)|=|g|*2`|g2. |πNovice 
 6372  programmers who calculate the standard deviation 
 6378  of some observations by using the textbook formula|'
 6386  {A9}|ε|≤s|4α=↓|4|↔K{H12}|↔a{H10}n|4|↔k|uc|)1|¬Ek|¬En|)x|=|βk
 6386  |g2α_↓|↔a|↔k|uc|)1|¬Ek|¬En|)x|βk|↔s|g2{H12}|↔s{H10}|↔hn(n|4α
 6386  _↓|41)|J!(14)|;{A9}*?*?*?|πoften _nd them selves 
 6391  taking the square root of a negative number*3 
 6399  A much better way to calculate means and standard 
 6408  deviations with ⊗oating-point arithmetic is to 
 6414  use the recurrence formulas|'{A9}|h|εS|β1|4|∂α=↓|40,!!S|βk|4
 6418  |∂α=↓|4S|βk|βα_↓|β1|4|↔N|4(x|βk|4|↔N|4M|βk|βα_↓|β1)|4|↔N|4(x
 6418  |βk|4|↔N|4M|βk),|E|n|;| M|β1|4|Lα=↓|4x|β1,| M|βk|4|Lα=↓|4M|β
 6419  k|βα_↓|β1|4|↔V|4(x|βk|4|↔B|4M|βk|βα_↓|β1)|4|↔n|4k,|J!(15)>
 6420  *?{A4}| S|β1|4|Lα=↓|40,| S|βk|4|Lα=↓|4S|βk|βα_↓|β1|4|↔V|4(x|β
 6420  k|4|↔B|4M|βk|βα_↓|β1)|4|↔N|4(x|βk|4|↔B|4M|βk),|J!(16)>
 6421  {A9}|πfor |ε2|4|¬E|4k|4|¬E|4n, |πwhere |ε|≤s|4α=↓|4|¬H|v4S|β
 6424  n/(n|4α_↓|41)|). |π[Cf. B. P. Welford, |εTechnometrics 
 6430  |≡4 (1962), 419<420.] |πWith this method |εS|βn 
 6437  |πcan never be negative, and we avoid other serious 
 6446  problems encountered by the naive method of accumulating 
 6454  sums, as shown in exercise 16. (See exercise 
 6462  19 for a summation technique that provides an 
 6470  even!|4|4|4|4Although algebraic laws do not always 
 6476  hold exactly, we can often show that they aren't 
 6485  too far o= base. When |εb|ge|gα_↓|g1|4|¬E|4x|4|¬W|4b|ge 
 6491  |πwe have round|ε(x)|4α=↓|4x|4α+↓|4|≤r(x) |πwhere 
 6495  |¬G|ε|≤r(x)|¬G|4|¬E|4|f1|d32|)b|ge|gα_↓|gp; |πhence|'
 6497  {A9}|πround(|εx)|4α=↓|4x{H12}({H10}1|4α+↓|4|≤d(x){H12}){H10}
 6497  ,|J!(17)|;{A9}|πwhere the relative error is bounded 
 6504  independently of |εx:|'{A9}|¬G|≤d(x)|¬G|4|¬E|4|f1|d32|)b|g1|
 6507  gα_↓|gp.|J!(18)|;{A9}|πWe can use this inequality 
 6513  to estimate the relative error of normalized 
 6520  ⊗oating-point calculations in a simple way, since 
 6527  |εu|4|↔V|4v|4α=↓|4(u|4α+↓|4v){H12}({H10}1|4α+↓|4|≤d(u|4α+↓|4
 6527  v){H12}){H10}, etc.|'|π!|4|4|4|4As an example 
 6532  of typical error-estimation procedures, let us 
 6538  consider the associative law for multiplication. 
 6544  Exercise 3 shows that |ε(u|4|↔N|4v)|4|↔N|4w |πis 
 6550  not in general equal to |εu|4|↔N|4(v|4|↔N|4w); 
 6556  |πbut the situation in this case is much better 
 6565  than it was with respect to the associative law 
 6574  of addition (1) and the distributive law (12). 
 6582  Infact, we have|'{A9}|ε(u|4|↔N|4v)|4|↔N|4w|4|∂α=↓|4{H12}({H1
 6585  0}(uv)(1|4α+↓|4|≤d|β1){H12}){H10}|4|↔N|4w|4α=↓|4uvw(1|4α+↓|4
 6585  |≤d|β1)(1|4α+↓|4|≤d|β2),|;{A4}| u|4|↔N|4(v|4|↔N|4w)|4|Lα=↓|4
 6586  u|4|↔N|4{H12}({H10}(vw)(1|4α+↓|4|≤d|β3){H12}){H10}|4α=↓|4uvw
 6586  (1|4α+↓|4|≤d|β3)(1|4α+↓|4|≤d|β4)>{A9}|πfor some 
 6589  |ε|≤d|β1, |≤d|β2, |≤d|β3, |≤d|β4, |πprovided 
 6594  that no exponent under⊗ow or over⊗ow occurs, 
 6601  where |ε|¬G|≤d|βj|¬G|4|¬E|4|f1|d32|)b|g1|gα_↓|gp 
 6603  |πfor each |εj. |πHence|'{A9}|(|ε(u|4|↔N|4v)|4|↔N|4w|d2u|4|↔
 6607  N|4(v|4|↔N|4w)|)|4α=↓|4|((1|4α+↓|4|≤d|β1)(1|4α+↓|4|≤d|β2)|d2
 6607  (1|4α+↓|4|≤d|β3)(1|4α+↓|4|≤d|β4)|)|4α=↓|41|4α+↓|4|≤d,|;
 6608  {A9}|πwhere|'{A9}|ε|¬G|≤d|¬G|4|¬E|42b|g1|gα_↓|gp/(1|4α_↓|4|f
 6609  1|d32|)b|g1|gα_↓|gp)|g2.|J!(19)|;{A9}|π!|4|4|4|4The 
 6611  number |εb|g1|gα_↓|gp |πoccurs so often in such 
 6618  analyses, it has been given a special name, one 
 6627  |εulp (|πmeaning one ``unit in the last place'' 
 6635  of the fraction part). Floating-point operations 
 6641  are correct to within half an |εulp, |πand the 
 6650  calculation of |εuvw |πby two ⊗oating-point multiplications 
 6657  will be correct within about one |εulp |π(ignoring 
 6665  second-order terms); hence the associative law 
 6671  for multiplication holds to within about two 
 6678  |εulps |πof relative error.|'!|4|4|4|4We have 
 6684  shown that |ε(u|4|↔N|4v)|4|↔N|4w |πis |εapproximately 
 6689  equal to u|4|↔N|4(v|4|↔N|4w), |πexcept when exponent 
 6695  over⊗ow or under⊗ow is a problem. It is worth 
 6704  while studying this intuitive idea of being ``approximately 
 6712  equal'' in more detail; can we make such a statement 
 6722  more precise in a reasonable way?|'!|4|4|4|4A 
 6729  programmer using ⊗oating-point arithmetic almost 
 6734  never wants to test if |εu|4α=↓|4v |π(or at least 
 6743  he hardly ever should try to do so), because 
 6752  this is an extremely improbable occurrence. For 
 6759  example, if a recurrence relation|'{A9}|εx|βn|βα+↓|β1|4α=↓|4
 6764  f|1(x|βn)|;{A9}|πis being used, where the theory 
 6771  in some textbook says |εx|βn |πapproaches a limit 
 6779  as |εn|4|¬M|4|¬X, |πit is usually a mistake to 
 6787  wait until |εx|βn|βα+↓|β1|4α=↓|4x|βn |πfor some 
 6792  |εn, |πsince the sequence |εx|βn |πmight be periodic 
 6800  with a longer period due to the rounding of intermediate 
 6810  results. The proper procedure is to wait until 
 6818  |¬G|εx|βn|βα+↓|β1|4α_↓|4x|βn|¬G|4|¬W|4|≤d, |πfor 
 6820  some suitably chosen number |ε|≤d; |πbut since 
 6827  we don't necessarily know the order of magnitude 
 6835  of |εx|βn |πin advance, it is even better to 
 6844  wait until|'{A9}|ε|¬Gx|βn|βα+↓|β1|4α_↓|4x|βn|¬G|4|¬E|4|≤e|¬G
 6846  x|βn|¬G;|J!(20)|;{A8}|πnow |ε|≤e |πis a number 
 6852  that is much easier to select. This relation 
 6860  (20) is another way of saying that |εx|βn|βα+↓|β1 
 6868  |πand |εx|βn |πare approximately equal; and our 
 6875  discussion indicates that a relation of ``approximately 
 6882  equal'' would be more useful than the traditional 
 6890  rel{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W.)!f
folio 293 galley 9
 6891  . 293!ch. 4!g. 8a|'{A12}!|4|4|4|4In other words, 
 6898  the fact that strict equality of ⊗oating-point 
 6905  values is of little importance implies that we 
 6913  ought to have a new operation, |ε⊃oating-point 
 6920  comparison, |πwhich is intended to help assess 
 6927  the relative values of two ⊗oating-point quantities. 
 6934  The following de_nitions seem to be appropriate, 
 6941  for base |εb, |πexcess |εq, |π⊗oating-point numbers 
 6948  |εu|4α=↓|4(e|βu,|4f|βu) |πand |εv|4α=↓|4(e|βv,|4f|βv):|'
 6951  |h|εu|4|↔Z|4v!(|≤e)!if|4and|4only|4if!|¬Gv|4α_↓|4u|¬G|4|∂|≤E
 6951  |4|≤e|4|πmax(|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|E|n|;
 6952  {A9}u|4|"2|4v!(|≤e)!|πif|4and|4only|4if| |εv|4α_↓|4u|4|L|¬Q|
 6952  4|≤e|4|πmax(|εb|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(21)|;
 6953  {A4}u|4|¬Z|4v!(|≤e)!|πif|4and|4only|4if!|¬G|εv|4α_↓|4u|¬G|4|
 6953  ¬E|4|≤e|4|πmax(|εb|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(22)|
 6953  ;{A4}u|4|"2|4v!(|≤e)!|πif|4and|4only|4if!|9|1|εu|4α_↓|4v|4|¬
 6954  Q|4|≤e|4|πmax|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq);|J!(23)|;
 6955  {A4}|εu|4|¬V|4v!(|≤e)!|πif|4and|4only|4if!|ε|¬Gv|4α_↓|4u|¬G|
 6955  4|¬E|4|≤e|4|πmin|ε(b|ge|ru|gα_↓|gq,|4b|ge|rv|gα_↓|gq).|J!(24
 6955  )|;{A9}|πThese de_nitions apply to unnormalized 
 6961  values as well as to normalized ones. Note that 
 6970  exactly less than), |εu|4|¬Z|4v (|πapproximately 
 6975  equal to), or |εu|4|"2|4v |π(de_nitely greater 
 6981  than) must always hold for any given pair of 
 6990  values |εu |πand |εv. |πThe relation |εu|4|¬V|4v 
 6997  |πis somewhat stronger than |εu|4|¬Z|4v, |πand 
 7003  it misht be read ``|εu |πis essentially equal 
 7011  to |εv.'' |πAll of the relations are given in 
 7020  terms of a positive real number |ε|≤e |πwhich 
 7028  measures the degree of approximation being considered.|'
 7035  !|4|4|4|4One way to view the above de_nitions 
 7042  is to associate a set |εS(u)|4α=↓|4|¬Tx|4|¬G|4|¬Gx|4α_↓|4u|¬
 7047  G|4|¬E|4|≤eb|ge|ru|gα_↓|gq|¬Y |πwith each ⊗oating-point 
 7051  number |εu; S(u) |πrepresents a set of values 
 7059  near |εu |πbased on the exponent of |εu'|πs ⊗oating-point 
 7068  representation. In these terms, we have |εu|4|"2|4v 
 7075  |πif and only if |εS(u)|4|¬W|4v |πand |εu|4|¬W|4S(v); 
 7082  u|4|¬Z|4v |πif and only if |εu|4|¬A|4S(v) |πor 
 7089  |εv|4|¬A|4S(u); u|4|"2|4v |πif and only if |εu|4|¬Q|4S(v) 
 7096  |πand |εS(u)|4|¬Q|4v; u|4|¬V|4v |πif and only 
 7102  if |εu|4|¬A|4S(v) |πand |εv|4|¬A|4S(u). |π(Here 
 7107  we are assuming that the parameter |ε|≤e, |πwhich 
 7115  measures the degree of approximation, is a constant; 
 7123  a more complete notation would indicate the dependence 
 7131  of |εS(u) |πupon |ε|≤e.)|'|π!|4|4|4|4Here are 
 7137  some simple consequences of the above de_nitions:|'
 7144  {A9}if!!|εu|4|"2|4v!(|≤e)!!|πthen!!|εv|4|"2|4u!(|≤e);|;
 7145  {A4}|πif!!|εu|4|¬V|4v!(|≤e)!!|πthen!!|εu|4|¬Z|4v!(|≤e);|J!(2
 7145  6)|;{A4}u|4|¬V|4u!(|≤e);|;{A4}|πif!!|εu|4|"2|4v!(|≤e)!!|πthe
 7147  n!!|εu|4|¬W|4v!|9|1|1|1|4|1|1|9|J!(28)|;{A6}|πif!!|εu|4|"2|4
 7148  v!(|≤e|β1)!|πand!|ε|≤e|β1|4|¬R|4|≤e|β2!!|πthen!!|εu|4|"2|4v!
 7148  (|≤e|β2);|J!(29)|;{A4}|πif!!|εu|4|¬Z|4v!(|≤e|β1)!|πand!|ε|≤e
 7149  |β1|4|¬E|4|≤e|β2!!|πthen!!|εu|4|¬Z|4v!(|≤e|β2);|J!(30)|;
 7150  {A4}|πif!!|εu|4|¬V|4v!(|≤e|β1)!|πand!|ε|≤e|β1|4|¬E|4|≤e|β2!!
 7150  |πthen!!|εu|4|¬V|4v!(|≤e|β2);|J!(31)|;{A4}|πif!!|εu|4|"2|4v!
 7151  (|≤e|β1)!|πand!|εv|4|"2|4w!(|≤e|β2)!!|πthen!!|εu|4|"2|4w!(|≤
 7151  e|β1|4α+↓|4|≤e|β2);|J!(32)|;{A4}|πif!!|εu|4|¬V|4v!(|≤e|β1)!|
 7152  πand!|εv|4|¬V|4w!(|≤e|β2)!!|πthen!!|εu|4|¬Z|4w!(|≤e|β1|4α+↓|
 7152  4|≤e|β2).|J!(33)|;{A9}|πMoreover, we can prove 
 7157  without di∃culty that|'{A9}|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬
 7160  G!|πand!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!!|εu|4
 7160  |¬V|4!(|≤e);|J!(34)|;{A4}|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬G!|π
 7161  and!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!!|εu|4|¬V|
 7161  4v!(|≤e);|J!(34)|;{A4}|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gu|¬G!|4|1|
 7162  1|πor!|4|1|1|1|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4|≤e|¬Gv|¬G!!|πimplies!
 7162  !|εu|4|¬Z|4v!(|≤e);|J!(35)|;{A9}|πand conversely, 
 7165  for |εnormalized |π⊗oating-point numbers |εu 
 7170  |πand |εv, |πwhen |ε|≤e|4|¬W|41,|'{A9}u|4|¬V|4v!(|≤e)!!|πimp
 7174  lies!!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gu|¬G!|πand!|ε|¬Gu|4α_↓|
 7174  4v|¬G|4|¬E|4b|≤e|¬Gv|¬G;|J!(36)|;{A4}u|4|¬Z|4v!(|≤e)!!|πimpl
 7175  ies!!|ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gu|¬G!|4|1|1|πor!|1|1|1|4
 7175  |ε|¬Gu|4α_↓|4v|¬G|4|¬E|4b|≤e|¬Gv|¬G.|J!(37)|;
 7176  {A9}|π|4|4|4|4!Let |ε|≤e|β0|4α=↓|4b|g1|gα_↓|gp 
 7178  |πone |εulp. |πFrom the derivation of (17) we 
 7186  have |ε|¬Gx|4α_↓|4|πround|ε(x)|¬G|4α=↓|4|¬G|≤r(x)|¬G|4|¬E|4|
 7187  f1|d32|)|≤e|β0|4|πmin(|ε|¬Gx|¬G, |¬G|πround|ε(x)|¬G), 
 7189  |πhence|'{A9}|εx|4|¬V|4|πround|ε(x)!(|f1|d32|)|≤e|β0);|J!(38
 7190  )|;{A9}|πand |εu|4|↔V|4v|4|¬V|4u|4α+↓|4v(|f1|d32|)|≤E|β0), 
 7193  |πetc. The approximate associative law for multiplication 
 7200  derived above can be recast as follows: We have|'
 7209  {A9}|¬G(|εu|4|↔N|4v)|4|↔N|4w|4α_↓|4u|4|↔N|4(v|4|↔N|4w)|¬G|4|
 7209  ¬E|4|(2|≤e|β0|d2(1|4α_↓|4|f1|d32|)|≤e|β0)|g2|)|4|¬Gu|4|↔N|4(
 7209  v|4|↔N|4w)|¬G|;{A9}|πby (19), and the same inequality 
 7216  is valid with (|εu|4|↔N|4v)|4|↔N|4w |πand |εu|4|↔N|4(v|4|↔N|
 7221  4w) |πinterchanged. Hence by (34),|'{A9}|ε(u|4|↔N|4v)|4|↔N|4
 7226  w|4|¬V|4u|4|↔N|4(v|4|↔N|4w)!(|≤e)|J!(39)|;{A9}|πwhenever 
 7228  |ε|≤e|4|¬R|42|≤e|β0/(1|4α_↓|4|f1|d32|)|≤e|β0)|g2. 
 7229  |πFor example, if |εb|4α=↓|410 |πand |εp|4α=↓|48 
 7235  |πwe may take |ε|≤e|4α=↓|40.00000021.|'|H*?{U0}{H10L12M29}583
folio 301 galley 10
 7239  20#Computer Programming!(Knuth/A.-W.)!f. 301!ch. 
 7242  4!g. 10a|'{A12}!|4|4|4|4The relations |"2, |¬Z, 
 7248  |"2, and |¬V are useful within numerical algorithms, 
 7256  and it is therefore a good idea to provide routines 
 7266  for comparing ⊗oating-point numbers as well as 
 7273  for doing arithmetic on them.|'!|4|4|4|4Let us 
 7280  now shift our attention back to the question 
 7288  of _nding |εexact |πrelations which are satis_ed 
 7295  by the ⊗oating-point operations. It is interesting 
 7302  to note that ⊗oating-point addition and subtraction 
 7309  are not completely intractable from an axiomatic 
 7316  standpoint, since they do satisfy the nontrivial 
 7323  identities stated in the following theorems:|'
 7329  {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡A|≡.|9|4|εLet u 
 7332  and v be normalized ⊃oating point numbers. Then|'
 7340  {A9}{H12}({H10}(u|4|↔V|4v)|4|↔B|4u{H12}){H10}|4α+↓|4{H12}({H
 7340  10}(u|4|↔V|4v)|4|↔B|4{H12}({H10}(u|4|↔V|4v)|4|↔B|4u){H12}){H
 7340  10}|4α=↓|4u|4|↔V|4v,|J!(40)|;{A9}{H10L12M29}provided 
 7342  that no exponent over⊃ow or under⊃ow occurs.|'
 7349  {A12}|πThis rather cumbersome looking identity 
 7354  can be rewritten in a simpler manner: Let|'{A9}|h|εu|¬C|4|∂α
 7362  =↓|4(u|4|↔V|4v)|4|↔B|4v|¬S,!!v|¬C|4|∂α=↓|4(u|4|↔V|4v)|4|↔B|4
 7362  u|¬S.|E|n|;{A8}(41)|E|?{B8}| u|¬S|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4
 7364  v,| v|¬S|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4u;>{A4}| u|¬C|4|Lα=↓|4(u|
 7365  4|↔V|4v)|4|↔B|4v|¬S,| v|¬C|4|Lα=↓|4(u|4|↔V|4v)|4|↔B|4u|¬S.>
 7366  {A9}|πIntuitively, |εu|¬S |πand |εu|¬C |πshould 
 7371  be approximations to |εu, |πand |εv|¬S |πand 
 7378  |εv|¬C |πshould be approximations to |εv. |πTheorem 
 7385  A tells us that|'{A9}|εu|4|↔V|4v|4α=↓|4u|¬S|4α+↓|4v|¬C|4α=↓|
 7389  4u|¬C|4α+↓|4v|¬S.|J!(42)>{A9}|πThis is a stronger 
 7394  statement than the identity|'{A9}|εu|4|↔V|4v|4α=↓|4u|¬S|4|↔V
 7398  |4v|¬C|4α=↓|4u|¬C|4|↔V|4v|¬S,|J!(43)|;{A9}|πwhich 
 7400  follows by rounding (42).|'{A12}|εProof.|9|4|πLet 
 7405  us say that |εt |πis a |εtail |πof |εx |πmodulo 
 7415  |εb|ge |πif|'{A9}|εt|4α=↓|4x(|πmodulo|4|εb|ge),!!|¬Gt|¬G|4|¬
 7417  E|4|f1|d32|)b|ge;|J!(44)|;{A9}|πthus, |εx|4α_↓|4|πround|ε(x)
 7419   |πis always a tail of |εx. |πThe proof of Theorem 
 7430  A rests largely on the following simple fact 
 7438  proved in exercise 11:|'{A12}|≡L|≡e|≡m|≡m|≡a 
 7443  |≡T.|4|9|εIf t is a tail of the ⊃oating-point 
 7451  number x then x|4|↔B|4t|4α=↓|4x|4α_↓|4t.|'{A12}|π!|4|4|4|4Le
 7455  t |εw|4α=↓|4u|4|↔V|4v. |πTheorem A holds trivially 
 7461  when |εw|4α=↓|40. |πBy multiplying all variables 
 7467  by a suitable power of |εb, |πwe may assume without 
 7477  loss of generality that |εe|βw|4α=↓|4p. |πThen 
 7483  |εu|4α+↓|4v|4α=↓|4w|4α+↓|4r, |πwhere |εr |πis 
 7487  a tail of |εu|4α+↓|4v |πmodulo 1. Furthermore 
 7494  |εu|¬S|4α=↓|4|πround(|εw|4α_↓|4v)|4α=↓|4|πround(|εu|4α_↓|4r)
 7494  |4α=↓|4u|4α_↓|4r|4α_↓|4t, |πwhere |εt |πis a 
 7499  tail of |εu|4α_↓|4r |πmodulo |εb|ge |πand |εe|4α=↓|4e|βu|4|4
 7505  α_↓|4p.|'|π!|4|4|4|4If |εe|4|¬E|40 |πthen |εt|4!|4u|4α_↓|4r|
 7509  4!|4|→α_↓v |π(modulo |εb|ge), |πhence |εt |πis 
 7515  a tail of |ε|→α_↓v |πand |εv|¬C|4α=↓|4|πround(|εw|4α_↓|4u|¬S
 7520  )|4α=↓|4|πround(|εv|4α+↓|4t)|4α=↓|4v|4α+↓|4t; 
 7521  |πthis proves (40). If |εe|4|¬Q|40 |πthen |ε|¬Gu|4α_↓|4r|¬G|
 7527  4|¬R|4b|gp|4α_↓|4|f1|d32|), |πand since |ε|¬Gr|¬G|4|¬E|4|f1|
 7530  d32|) |πwe have |ε|¬Gu|¬G|4|¬R|4b|gp|4α_↓|41. 
 7534  |πIt follows that |εr |πis a tail of |εv |πmodulo 
 7544  1. If |εu|¬S|4α=↓|4u, |πwe have |εv|¬C|4α=↓|4|πround(|εw|4α_
 7549  ↓|4u)|4α=↓|4|πround(|εv|4α_↓|4r)|4α=↓|4v|4α_↓|4r. 
 7550  |πOtherwise the relation round(|εu|4α_↓|4r)|4|=|↔6α=↓|4u 
 7554  |πimplies that |ε|¬Gu|¬G|4α=↓|4b|gp|4α_↓|41, 
 7557  |¬Gr|¬G|4α=↓|4|f1|d32|), |¬Gu|¬S|¬G|4α=↓|4b|gp; 
 7559  |πwe have |εv|¬C|4α=↓|4|πround|ε(w|4α_↓|4u|¬S)|4α=↓|4|πround
 7561  (|εv|4α+↓|4r)|4α=↓|4v|4α+↓|4r |πbecause |εr |πis 
 7565  also a tail of |ε|→α_↓v |πin this case|'{A12}!|4|4|4|4Theore
 7573  m A exhibits a regularity property of ⊗oating-point 
 7581  addition, but it doesn't seem to be an especially 
 7590  useful result. The following identity is more 
 7597  signi_cant:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡B.|9|4|εUnder 
 7600  the hypotheses of Theorem A and (41),|'{A9}u|4α+↓|4v|4α=↓|4(
 7607  u|4|↔V|4v)|4α+↓|4{H12}({H10}(u|4|↔B|4u|¬S)|4|↔V|4(v|4|↔B|4v|
 7607  ¬C){H12}){H10}.|J!(45)|;{A12}Proof.|9|4|πIn fact, 
 7610  we can show that |εu|4|↔B|4u|¬S|4α=↓|4u|4α_↓|4u|¬S, 
 7615  v|4|↔B|4v|¬C|4α=↓|4v|4α_↓|4v|¬C, |πand |ε(u|4α_↓|4u|¬S)|4|↔V
 7617  |4(v|4α_↓|4v|¬C)|4α=↓|4(u|4α_↓|4u|¬S)|4α+↓|4(v|4α_↓|4v|¬C), 
 7618  |πhence (45) will follow from Theorem A. Using 
 7626  the notation of the preceding proof, these relations 
 7634  are respectively equivalent to|'{A9}round(|εt|4α+↓|4r)|4α=↓|
 7638  4t|4α+↓|4r,!!|πround|ε(t)|4α=↓|4t,!!|πround|ε(r)|4α=↓|4r.|J!
 7638  (46)|;{A9}|πExercise 12 establishes the theorem 
 7644  in the special case |ε|¬Ge|βu|4α_↓|4e|βv|¬G|4|¬R|4p. 
 7649  |πOtherwise |εu|4α+↓|4v |πhas at most 2|εp |πsigni_cant 
 7656  digits and it is easy to see that round(|εr)|4α=↓|4r. 
 7665  |πIf now |εe|4|¬Q|40, |πthe proof of Theorem 
 7672  A shows that |εt|4α=↓|4|→α_↓r |πor |εt|4α=↓|4r|4α=↓|4|→|¬N|f
 7677  1|d32|). |πIf |εe|4|¬E|40 |πwe have |εt|4α+↓|4r|4!|4u 
 7683  |πand |εt|4!|4|→α_↓v (|πmodulo|4|εb|ge); |πthis 
 7687  is enough to prove that |εt|4α+↓|4r |πand |εr 
 7695  |πround to themselves, provided that |εe|βu|4|¬R|4e 
 7701  |πand |εe|βv|4|¬R|4e. |πBut either |εe|βu|4|¬W|40 
 7706  |πor |εe|βv|4|¬W|40 |πwould contradict our hypothesis 
 7712  that |ε|¬Ge|βu|4α_↓|4e|βv|¬G|4|¬W|4p, |πsince 
 7715  |εe|βw|4α=↓|4p.|'|π!|4|4|4|4Theorem B gives |εan 
 7720  explicit formula for the di∩erence |πbetween 
 7726  |εu|4α+↓|4v |πand |εu|4|↔V|4v, |πin terms of 
 7732  quantities that can be calculated directly using 
 7739  _ve operations of ⊗oating-point arithmetic. If 
 7745  the radix |εb |πis 2 or 3, we can improve on 
 7756  this result, obtaining the exact value of the 
 7764  correction term with only two ⊗oating-point operations 
 7771  and one (_xed-point) comparison of absolute values:|'
 7778  {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡C|≡.|9|4|εIf b|4|¬E|43 
 7781  and |¬Gu|¬G|4|¬R|4|¬Gv|¬G then|'{A9}u|4α+↓|4v|4α=↓|4(u|4|↔V|
 7784  4v)|4α+↓|4{H12}({H10}u|4|↔B|4(u|4|↔V|4v){H12}){H10}|4|↔V|4v.
 7784  |J!(47)|;{A12}Proof.|9|4|πFollowing the conventions 
 7788  of preceding proofs again, we wish to show that 
 7797  |εv|4|↔B|4v|¬S|4α=↓|4r. |πIt su∃ces to show that 
 7803  |εv|¬S|4α=↓|4w|4α_↓|4u, |πbecause (46) will then 
 7808  yield |εv|4|↔B|4v|¬S|4α=↓|4|πround(|εv|4α_↓|4v|¬S)|4α=↓|4|πr
 7809  ound(|εu|4α+↓|4v|4α_↓|4w)|4α=↓|4|πround|ε(r)|4α=↓|4r.|'
 7810  !|4|4|4|4|πWe shall in fact prove (47) whenever 
 7817  |εb|4|¬E|43 |πand |εe|βu|4|¬R|4e|βr. |πIf |εe|βu|4|¬R|4p 
 7822  |πthen |εr |πis a tail of |εv |πmodulo 1, hence 
 7832  |εv|¬S|4α=↓|4w|4|↔B|4u|4α=↓|4v|4|↔B|4r|4α=↓|4v|4α_↓|4r|4α=↓|
 7832  4w|4α_↓|4u |πas desired. If |εe|βu|4|¬W|4p |πthen 
 7838  we must have |εe|βu|4α=↓|4p|4α_↓|41, |πand |εw|4α_↓|4u 
 7844  |πis a multiple of |εb|gα_↓|g1; |πit will therefore 
 7852  round to itself if its magnitude is less than 
 7861  |εb|gp|gα_↓|g1|4α+↓|4b|gα_↓|g1. |πSince |εb|4|¬E|43, 
 7864  |πwe have indeed |ε|¬Gw|4α_↓|4u|¬G|4|¬E|4|¬Gw|4α_↓|4u|4v|¬G|
 7867  4α+↓|4|¬Gv|¬G|4|¬E|4|f1|d32|)|4α+↓|4(b|gp|gα_↓|g1|4α_↓|4b|gα
 7867  _↓|g1)|4|¬W|4b|gp|gα_↓|g1|4α+↓|4b|gα_↓|g1.|'{A12}|π!|4|4|4|4
 7868  The proofs of Theorems A, B, and C do not rely 
 7879  on the precise de_nitions of round(|εx) |πin 
 7886  the ambiguous cases when |εx |πis exactly midway 
 7894  between consecutive ⊗oating-point numbers; any 
 7899  way of resolving the ambiguity will su∃ce for 
 7907  the validity of everything we have proved so 
 7915  far. No rounding rule can be best for every application; 
 7925  for example, we generally want a special rule 
 7933  when computing our income tax. But for most numerical 
 7942  calculations the best policy appears to be the 
 7950  rounding scheme speci_ed in Algorithm 4.2.1 N, 
 7957  which insists that the least signi_cant digit 
 7964  should always be made even (or always odd) when 
 7973  an ambiguous value is rounded. This is not a 
 7982  trivial technicality, of interest only to nit-pickerss; 
 7989  it is an important practical consideration, since 
 7996  the ambiguous case arises surprisingly often 
 8002  and a biased rounding rule produces signi_cantly 
 8009  poor results. For example, consider decimal arithmetic 
 8016  and assume that remainders of 5 are always rounded 
 8025  upwards. Then if |εu|4α=↓|41.0000000 |πand |εv|4α=↓|40.55555
 8030  555 |πwe have |εu|4|↔V|4v|4α=↓|41.5555556; |πand 
 8035  if we ⊗oating-subtract |εv |πfrom this result 
 8042  we get |εu|¬S|4α=↓|41.0000001. |πAdding and subtracting 
 8048  |εv |πfrom |εu|¬S |πgives 1.0000002, and the 
 8055  next time we will get 1.000003, etc.; the result 
 8064  keeps growing although we are adding and subtracting 
 8072  the same value*3|'!|4|4|4|4This phenomenon, called 
 8078  |εdrift, |πwill not occur when we use a stable 
 8087  rounding rule based on the parity of the least 
 8096  signi_cant digit. More precisely:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡
 8100  m |≡D|≡.|9|4|ε{H12}({H10}()u|4|↔V|4v)|4|↔B|4v)|4|↔V|4v{H12})
 8101  {H10}|4|↔B|4v|4α=↓|4(u|4|↔V|4v)|4|↔B|4v.|'{A12}|πFor 
 8103  example, if |εu|4α=↓|41.2345679 |πand |εv|4α=↓|4α_↓0.2345678
 8107  5, |πwe _nd |εu|4|↔V|4v|4α=↓|41.0000000, (u|4|↔V|4v)|4|↔B|4v
 8111  |4α=↓|41.2345678, {H12}({H10}(u|4|↔V|4v)|4|↔B|4v{H12}){H10}|
 8112  4|↔V|4v|4α=↓|40.99999995, {H12}({H10}((u|4|↔V|4v)|4|↔B|4v)|4
 8113  |↔V|4v{H12}){H10}|4|↔B|4v|4α=↓|41.2345678. |πThe 
 8115  proof for general |εu |πand |εv |πseems to require 
 8124  a case analysis even more detailed than that 
 8132  in the above theorems; see the references at 
 8140  the end of this section.|'|H*?*?*?*?{U0}{H10L12M29}58320#Compute
folio 300 galley 11
 8145  r Programming!(Knuth/A.-W.)!f. 300!ch. 4!g. 11a|'
 8150  {A12}!|4|4|4|4Theorem D is valid both for ``round 
 8157  to even'' and ``round to odd''; how should we 
 8166  choose between these possibilities? When the 
 8172  radix |εb |πis odd, ambiguous cases never arise 
 8180  in _nite precision except during ⊗oating-point 
 8186  division, and the rounding in such cases is comparatively 
 8195  unimportant. Experience shows that the best rule 
 8202  for |εeven |πradices is, ``Round to even when 
 8210  |εb/2 |πis odd, round to odd when |εb/2 |πis 
 8219  even.'' The least signi_cant digit of a ⊗oating-point 
 8227  fraction occurs frequently as a remainder to 
 8234  be rounded o= in subsequent calculations, and 
 8241  this rule avoids generating the digit |εb/2 |πin 
 8249  the least signi_cant position whenever possible; 
 8255  its e=ect is to provide some memory of an ambiguous 
 8265  rounding so that subsequent rounding will tend 
 8272  to be unambiguous. For example, if we were to 
 8281  round to odd in the decimal system, repeated 
 8289  rounding of the number 2.44445 to one less place 
 8298  each time leads to the sequence 2.4445, 2.445, 
 8306  2.45, 2.5, 3; |πbut if we round to even, such 
 8316  situations do not occur. [Roy A. Keir, |εInf. 
 8324  Proc. Letters |≡3 (1975), 188<189.]|'|π!|4|4|4|4A 
 8330  reader who has checked some of the details of 
 8339  the above proofs will realize the immense simpli_cation 
 8347  that has been a=orded by the simple rule |εu|4|↔V|4v|4α=↓|4|
 8355  πround(|εu|4α+↓|4v). |πIf our ⊗oating-point addition 
 8360  routine would tail to give this result even in 
 8369  a few rare cases, the proofs would become enormously 
 8378  more complicated and perhaps they would even 
 8385  break down completely.|'!|4|4|4|4Theorem B faills 
 8391  if truncation arithmetic is used in place of 
 8399  rounding, i.e., if we let |εu|4|↔V|4v|4α=↓|4|πtrunc|ε(u|4α+↓
 8404  |4v) |πand |εu|4|↔B|4v|4α=↓|4|πtrunc(|εu|4α_↓|4v), 
 8407  |πwhere trunc(|εx) |πtakes all positive real 
 8413  |εx |πinto the largest ⊗oating-point number |→|¬E|εx. 
 8420  |πAn exception to Theorem B would then occur 
 8428  for cases such as (20, |→α+↓.10000001)|4|↔V|4 
 8434  (10, |→α_↓.10000001)|4α=↓|4(20,|4|→α+↓.10000000), 
 8436  when the di=erence between |εu|4α+↓|4v |πand 
 8442  |εu|4|↔V|4v |πcannot be expressed exactly as 
 8448  a ⊗oating-point number; and also for cases such 
 8456  as 12345678|4|↔V|4.012345678, when it can be.|'
 8462  !|4|4|4|4Many people feel that, since ⊗oating-point 
 8468  arithmetic is inexact by nature, there is no 
 8476  harm in making it just a little bit less exact 
 8486  in certain rather rare cases, if it is convenient 
 8495  to do so. This policy saves a few cents in the 
 8506  design of computer hardware, or a small percentage 
 8514  of the average running time of a subroutine. 
 8522  But the above discussion shows that such a policy 
 8531  is mistaken. We could save about _ve percent 
 8539  of the running time of the |¬f|¬a|¬d|¬d subroutine, 
 8547  Program 4.2.1A, and about 25 percent of its space, 
 8556  if we took the liberty of rounding incorrectly 
 8564  in a few cases, but we are much better o= leaving 
 8575  it as it is. The reason is not to glorify ``bit 
 8586  chasing'' fundamental issue is at stake here: 
 8593  |εNumerical subroutines should deliver results 
 8598  which satisfy simple, useful mathematical laws 
 8604  whenever possible. |πThe crucial formula |εu|4|↔V|4v|4α=↓|4|
 8609  πround(|εu|4α+↓|4v) |πis a ``regularity'' property 
 8614  which makes a great deal of di=erence between 
 8622  whether mathematical analysis of computational 
 8627  algorithms is worth doing or worth avoiding*3 
 8634  Without any underlying symmetry properties, the 
 8640  job of proving interesting results becomes extremely 
 8647  unpleasant. |εThe enjoyment of one's tools is 
 8654  an essential ingredient of successful work.|'
 8660  {A12}|π|≡B|≡.|9|≡U|≡n|≡n|≡o|≡r|≡m|≡a|≡l|≡i|≡z|≡e|≡d 
 8661  |≡⊗|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡o|≡i|≡n|≡t |**≡a|π|≡B|≡.|9|≡U|≡n
 8662  |≡n|≡o|≡r|≡m|≡a|≡l|≡i|≡z|≡e|≡d |≡⊗|≡o|≡a|≡t|≡i|≡n|≡g|≡-|≡p|≡
 8663  o|≡i|≡n|≡t |≡a|≡r|≡i|≡t|≡h|≡m|≡e|≡t|≡i|≡c|≡.|9|4The 
 8665  policy of normalizing all ⊗oating-point numbers 
 8671  may be construed in two ways: We may look on 
 8681  it favorably by saying it is an attempt to get 
 8691  the maximum possible accuracy obtainable with 
 8697  a given degree of precision, or we may consider 
 8706  it to be potentially dangerous since it tends 
 8714  to imply that the results are more accurate than 
 8723  they really are. When we normalize the result 
 8731  of (1, |→α+↓.31428571)|4|↔B|4(1, |→α+↓.31415927) 
 8735  to (|→α_↓2, |→α+↓.12644000), we are suppressing 
 8741  information about the possibly greater inaccuracy 
 8747  of the latter quantity. Such information would 
 8754  be retained if the answer were left as (1, |→α+↓.00012644).|
 8763  '!|4|4|4|4The input data to a problem is frequently 
 8772  not known as precisely as the ⊗oating-point representation 
 8780  allows. For example, the values of Avogadro's 
 8787  number and Planck's constant are not known to 
 8795  eight signi_cant digits, and its might be more 
 8803  appropriate to denote them, respectively, by|'
 8809  {A9}(27,|4|→α+↓.00060225)!!and!!(|→α_↓23,|4|→α+↓.00010545)|;
 8810  {A9}instead of by (24, |→α+↓.60225200) and (|→α_↓26, 
 8817  |→α+↓.10545000). It would be nice if we could 
 8825  give our input data for each problem in an unnormalized 
 8835  form which expresses how much precison is assumed, 
 8843  and if the output would indicate just how much 
 8852  precision is known in the naswer. Unfortunately, 
 8859  this is a terribly di∃cult proble, although the 
 8867  use of unnormalized arithmetic can help to give 
 8875  some indication. For example, we can say with 
 8883  a fair degree of certainty that the product of 
 8892  Avogadro's number by Planck's constant is (0,|4|→α+↓.0006350
 8898  7), and their sum is (27, |→α+↓.00060225). (The 
 8906  purpose of this example is not to suggest that 
 8915  any important physical signi_cance should be 
 8921  attached to the sum and product of these fundamental 
 8930  constants; the point is that it is possible to 
 8939  preserve a little of the information about precision 
 8947  in the result of calculation with imprecise quantities, 
 8955  when the original operands are independent of 
 8962  each other.)|'!|4|4|4|4The rules for unnormalized 
 8968  arithmetic are simply this: Let |εl|βu |πbe the 
 8976  number of leading zeros in the fraction part 
 8984  of |εu|4α=↓|4(e|βu, f|βu), |πso that |εl|βu |πis 
 8991  the largest integer |ε|→|¬Ep |πwith |ε|¬G|1f|βu|¬G|4|¬W|4b|g
 8996  α_↓|gl|ru. |πThen addition and subtraction are 
 9002  performed just as in Algorithm 4.2.1A, except 
 9009  that all scaling to the left is suppressed. Multiplication 
 9018  and division are performed as in Algorithn 4.2.1M, 
 9026  except that the answer is scaled right or left 
 9035  so that precisely max (|εl|βu, l|βv) |πleading 
 9042  zeros appear. Essentially the same rules have 
 9049  been used in manual calculation for may years.|'
 9057  !|4|4|4|4It follows that, for unnormalized computations,|'
 9063  {A9}|h|εe|βu|β|↔V|βv,|4e|βu|β|↔B|βv|∂α=↓|4e|βu|4α_↓|4e|βv|4α
 9063  +↓|4q|4α_↓|4l|βu|4α+↓|4l|βv|4α+↓|4|πmax(|εl|βu,|4l|βv)|4α+↓|
 9063  4(0|4|πor|41).|E|n|;| |εe|βu|β|↔V|βv,|4e|βu|β|↔B|βv|4|Lα=↓|4
 9064  |πmax(|εe|βu,|4e|βv)|4α+↓|4(0|4|πor|41)|J!(48)>
 9065  {A5}| |εe|βu|β|↔N|βv|4|Lα=↓|4e|βu|4α+↓|4e|βv|4α_↓|4q|4α_↓|4|
 9065  πmin(|εl|βu,|4l|βv)|4α_↓|4(0|4|πor|41)|J!(49)>
 9066  {A4}| |εe|βu|β|↔n|βv|4|Lα=↓|4e|βu|4α_↓|4e|βv|4α+↓|4q|4α_↓|4l
 9066  |βu|4α+↓|4l|βv|4α+↓|4|πmax(|εl|βu,|4l|βv)|4α+↓|4(0|4|πor|41)
 9066  .|J!(50)>{A9}When the result of a calculation 
 9073  is zero, an unnormalized zero (often called an 
 9081  ``order of magnitude zero'') is given as the 
 9089  answer; this indicates that the answer may not 
 9097  truly be zero, we just don't known any of its 
 9107  signi_cant digits*3|'!|4|4|4|4Error analysis takes 
 9112  a somewhat di=erent form with unnormalized ⊗oating-point 
 9119  arithmetic. Let us de_ne|'{A9}|ε|≤d|βu|4α=↓|4|f1|d32|)b|ure|
 9123  βuα_↓qα_↓p|)|)!!|πif!!|εu|4α=↓|4(e|βu,|4f|βu).|J!(51)|;
 9124  {A9}|πThis quantity depends on the representation 
 9130  of |εu, |πnot just on the value |εb|ge|ru|gα_↓|gqf|βu. 
 9138  |πOur rounding rule tells us that|'*1{A9}|ε|¬Gu|4|↔V|4v|4α_↓|
 9144  4(u|4α+↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔V|βv,!!|¬Gu|4|↔B|4v|4α_↓|4(u
 9144  |4α_↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔B|βv,|;{A5}|¬Gu|4|↔N|4v|4α_↓|4(
 9145  u|4α⊗↓|4v)|¬G|4|¬E|4|≤d|βu|β|↔N|βv,!!|¬Gu|4|↔n|4v|4α_↓|4(|4|
 9145  4/|4|1|4v)|¬G|4|¬E|4|≤d|βu|β|↔n|βv.|;{A9}|πThese 
 9147  inequalities apply to normalized as well as unnormalized 
 9155  arithmetic; the main di=erence between the two 
 9162  types of error analysis is the de_nition of the 
 9171  exponent of the result of each operation (Eqs. 
 9179  (48) to (50){H12}){H10}.|'!|4|4|4We have remarked 
 9185  that the relations |"2, |¬Z, |"2, and |¬V de_ned 
 9194  earlier in this section are valid and meaningful 
 9202  for unnormalized numbers as well as for normalized 
 9210  numbers. As an example of the use of these relations, 
 9220  let us prove an approximate associative law for 
 9228  unnormalized addition, analogous to (39):|'{A9}|ε(u|4|↔V|4v)
 9233  |4|↔V|4w|4|¬V|4u|4|↔V|4(v|4|↔V|4w)!!(|≤e),|J!(52)|;
 9234  {A9}|πfor suitable |ε|≤e. |πWe have|'{A9}|ε|¬G(u|4|↔V|4v)|4|
 9239  ↔V|4w|4α_↓|4(u|4α+↓|4v|4α+↓|4w)|¬G|4|∂|¬E|4|¬G(u|4|↔V|4v)|4|
 9239  ↔V|4w|4α_↓|4{H12}({H10}(u|4|↔V|4v)|4α+↓|4w{H12}){H10}|¬G|;
 9240  {A4}|L!|4α+↓|4|¬Gu|4|↔V|4v|4α_↓|4(u|4α+↓|4v)|¬G|'
 9241  {A4}|L|¬E|4|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w|)α+↓|4|≤d|βu|1|1|β|↔V
 9241  |1|1|βv>{A4}|L|¬E|42|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w|).>
 9243  {A9}|πA similar formula holds for |ε|¬Gu|4|↔V|4(v|4|↔V|4w)|4
 9248  α_↓|4(u|4α+↓|4v|4α+↓|4w)|¬G. |πNow since |εe|ur|)(u|4|↔v|4v)
 9251  |4|↔V|4w|)|4α=↓|4|πmax|ε(e|βu,|4e|βv,|4e|βw)|4α+↓|4(0,|41,|4
 9251  |πor|42), we have |ε|≤d|ur|)(u|4|↔V|4v)|4|↔V|4w)|)|4|¬E|4b|g
 9254  2|≤d|ur|)u|4|↔V|4(v|4|↔V|4w)|). |πTherefore we 
 9257  _nd that (52) is valid when |ε|≤e|4|¬R|42b|g2|gα_↓|gp; 
 9264  |πunnormalized addition is not as erratic as 
 9271  normalized addition with respect to the associative 
 9278  law.|'|Hβ*?*?*?{U0}{H10L12M29}58320#Computer Programming!(Knuth
folio 306 galley 12
 9280  /A.-W.)!ch. 4!f. 306!g. 12a|'{A12}!|4|4|4|4It 
 9285  should be emphasized that unnormalized arithmetic 
 9291  is by no means a panacea; there are examples 
 9300  where it indicates greater accuracy than is present 
 9308  (e.g., adding up a great many small quantities 
 9316  of about the same magnitude, or _nding the |εn|πth 
 9325  power of a number for large |εn), |πand there 
 9334  are many more examples when it indicataes poor 
 9342  accuracy while normalize arithmetic actually 
 9347  does produce good results. There is an important 
 9355  reason why no straightforward one-operation-at-a-time 
 9360  method of error analysis is really reliable, 
 9367  namely the fact that operands are usually not 
 9375  independent of each other. This means that errors 
 9383  tend to cancel or reinforce each other in strange 
 9392  ways. For example, suppose that |εx |πis approximately 
 9400  |f1|d32|), and suppose that we have an approximation 
 9408  |εy|4α=↓|4x|4α+↓|4|≤d |πwith absolute error |ε|≤d. 
 9413  |πIf we now wish to compute |εx(1|4α_↓|4x), |πwe 
 9421  can form |εy(1|4α_↓|4y); |πif |εx|4α=↓|4|f1|d32|)|4α+↓|4|≤e 
 9426  |πwe _nd |εy(1|4α_↓|4y)|4α=↓|4x(1|4α_↓|4x)|4α_↓|42|≤e|≤d|4α_
 9428  ↓|4|≤d|g2: |πThe error has decreased by a factor 
 9436  of 2|ε|≤e|4α+↓|4|≤d. |πThis is just one case 
 9443  where multiplication of imprecise quantities 
 9448  can lead to a quite accurate result when the 
 9457  operands are not independent of each other. A 
 9465  more obvious example is the computation of |εx|4|↔B|4x, 
 9473  |πwhich can be obtained with perfect accuracy 
 9480  regardless of how bad an approximation to |εx 
 9488  |πwe begin with.|'!|4|4|4|4The extra information 
 9494  which unnormalized arithmetic gives us can often 
 9501  be more important than the information it destroys 
 9509  during an extended calculation, but (as usual) 
 9516  we must use it with care. Examples of the proper 
 9526  use of unnormalized arithmetic are discussed 
 9532  by R. L. Ashenhurst and N. Metropolis, in |εComputers 
 9541  and Computing, AMM |πSlaught Memorial Papers, 
 9547  |≡1|≡0 (February, 1965), 47<59; and by R. L. 
 9555  Ashenhurst, in |εError in Digital Computation, 
 9561  |πvol. II, ed. by L. B. Rall (New York: Wiley, 
 9571  1965), 3<37. Appropriate methods for computing 
 9577  standard mathematical functions with both input 
 9583  and output in unnormalized form are given by 
 9591  R. L. Ashenhurst in |εJACM |≡1|≡1 |π(1964), 168<187.|'
 9599  |≡C|≡.|9|≡I|≡n|≡t|≡e|≡r|≡n|≡a|≡l |≡a|≡r|≡i|≡t|≡h|≡m|≡e|≡t|≡i
 9600  |≡c|≡.|9|4Another approach to the problem of 
 9606  error determination is the so-called ``interval'' 
 9612  or ``range'' arithmetic, in which upper and lower 
 9620  bounds on each number are maintained during the 
 9628  calculations. Thus, for example, if we know that 
 9636  |εu|β0|4|¬E|4u|4|¬E|4u|β1 |πand |εv|β0|4|¬E|4v|4|¬E|4v|β1, 
 9639  we represent this by the interval notation |εu|4α=↓|4[u|β0,|
 9646  4u|β1], v|4α=↓|4[v|β0,|4v|β1]. |πThe sum |εu|4|↔V|4v 
 9651  |πis [|εu|β0|4!|4v|β0, u|β1|4!|4v|β1], |πwhere 
 9655  ! denotes ``lower ⊗oating-point addition,'' the 
 9661  greatest representable number less than or equal 
 9668  to the true sum, and ! is de_ned similarly (see 
 9678  exercise 4.2.1<13). Furthermore |εu|4|↔B|4v|4α=↓|4[u|β0|4!|4
 9681  v|β1, u|β1|4!|4v|β0]; |πand if |εu|β0 |πand |εv|β0 
 9688  |πare positive, we have |εu|4|↔N|4v|4α=↓|4[u|β0|4!|4v|β0, 
 9693  u|β1|4!|4v|β1], u|4|↔n|4v|4α=↓|4[u|β0|4!|4v|β1, 
 9695  u|β1|4!|4v|β0]. |πFor example, we might represent 
 9701  Avogadro's number and Planck's constant as|'{A9}|εN|4α=↓|4{H
 9707  12}[{H10}(24,|4|→α+↓.60222400),|4(24,|4|→α+↓.60228000){H12}]
 9707  {H10}, h|4α=↓|4{H12}[{H10}(|→α_↓26,|4|→α+↓.10544300),|4(|→α_
 9708  ↓26,|4|→α+↓.10545700){H12}]{H10};|;{A9}|πtheir 
 9710  sum and product would then turn out ot be|'{A9}|εN|4|↔V|4h|4
 9719  α=↓|4{H12}[{H10}(24,|4|→α+↓.60222400),|4(24,|4|→α+↓.60228001
 9719  ){H12}]{H10},|4N|4|↔N|4h|4α=↓|4{H12}[{H10}(|→α_↓3,|4|→α+↓.63
 9719  500305),|4(|→α_↓3,|4|→α+↓.63514642){H12}]{H10}.|;
 9720  {A9}{A9}|π!|4|4|4|4If we try to divide by [|εv|β0,|4v|β1] 
 9727  |πwhen |εv|β0|4|¬W|40|4|¬W|4v|β1, |πthere is 
 9731  a possibility of division by zero. Since the 
 9739  philosophy underlying interval arithmetic is 
 9744  to provide rigorous error estimates, a divide-by-zero 
 9751  error should be signalled in this case. However, 
 9759  over⊗ow and under⊗ow need not be treated as errors 
 9768  in interval arithmetic, if special conventions 
 9774  are introduced as discussed in exercise 24.|'
 9781  !|4|4|4|4Interval arithmetic takes only about 
 9786  twice as long as ordinary arithmetic, and it 
 9794  provides truly reliable error estimates. Considering 
 9800  the di∃culty of mathematical error analyses, 
 9806  this is indeed a small price to pay. Since the 
 9816  intermediate values in a calculation often depend 
 9823  on each other, as explained above, the _anal 
 9831  estimates obtained with interval arithmetic will 
 9837  tend to be pessimistic; and iterative numerical 
 9844  methods often have to be redesigned if we want 
 9853  to deal with intervals. The prospects for e=ective 
 9861  use of interval arithmetic look very good, however, 
 9869  and e=orts should be made to increase it availability.|'
 9878  {A12}|≡D|≡. |≡H|≡i|≡s|≡t|≡o|≡r|≡y |≡a|≡n|≡d |≡b|≡i|≡b|≡l|≡i|
 9881  ≡o|≡g|≡r|≡a|≡p|≡h|≡y|≡.|9|4Jules Tannery's classic 
 9884  treatise on decimal calculations, |εle|=&c0ns 
 9889  d'Arithm|=1etique |π(Paris: Colin, 1894), stated 
 9894  that positive numbers should be rounded upeards 
 9901  if the _rst discarded digit is 5 or more; since 
 9911  exactly half of the decimal digits are 5 or more, 
 9921  he felt that this rule rounds exactly half the 
 9930  time and produces compensating errors. The idea 
 9937  of ``round to even'' in the ambiguous cases seems 
 9946  to have been mentioned _rst by James B. Scarborough 
 9955  in the _rst edition of his pioneering book |εNumerical 
 9964  Mathematical Analysis |π(Baltimore: Uohns Hopkins 
 9969  Press, 1930), p. 2; in the second (1950) edition 
 9978  he ampli_ed his earlier remarks, stating that 
 9985  ``It should be obvious to any thinking person 
 9993  that when a 5 is cut o=, the preceding digit 
10003  should be increased by 1 in only |εhalf |πthe 
10012  cases,'' and he recommended round-to-even in 
10018  order to achieve this.|'!|4|4|4|4The _rst analysis 
10025  of ⊗oating-point arithmetic was given by F. L. 
10033  Bauer and K. Samelson, |εZeitschrift f|=4ur angewandte 
10040  Math. und Physik |≡5 (1953), 312<316. |πThe next 
10048  publication was not until over _ve years later: 
10056  J. W. Carr III, |εCACM |≡2, 5 |π(May, 1959), 
10065  10<15. See also P. C. Fischer, |εProc. ACM |*/|↔O|↔L|\th 
10074  Nat. Meeting (|πUrbana, Illinois, 1958), paper 
10080  39. The book |εRounding Errors in Algebraic Processes 
10088  (|πEnglewood Cli=s: Prentice-Hall, 1963), by 
10093  J. H. Wilkinson, shows how to apply error analysis 
10102  of the individual arithmetic operations to the 
10109  error analysis of large-scale problems; see also 
10116  his treatise on |εThe Algebraic Eigenvalue Problem 
10123  (|πOxford: Clarendon Press, 1965).|'!|4|4|4|4The 
10128  relations |"2, |¬Z, |"2, |¬V introduced in this 
10136  section are similar to ideas published by A. 
10144  van Wijngaarden in |εBIT |≡6 (1966), 66<81. |πTheorems 
10152  A and B above were inspired by some related work 
10162  of Ole M|=|↔6oller, |εBIT |≡5 (1965), 37<50, 
10169  251<255. |πTheorem C is due to T. J. Dekker, 
10178  |εNumer. Math. |≡1|≡8 (1971), 224<242; |πextensions 
10184  and re_nements of all three theorems have been 
10192  published by S. Linnainmaa, |εBIT |π|≡1|≡4 (1974), 
10199  167<202. W. M. Kahan introduced Theorem D in 
10207  some unpublished notes; for a complete proof 
10214  and further commentary see J. F. Reiser and D. 
10223  E. Knuth, |εInf. Proc. Letters |≡3 (1975), 84<87, 
10231  164.|'|π!|4|4|4|4Unnormalized ⊗oating-point arithmetic 
10235  was recommended by F. L. Bauer and K. Samelson 
10244  in the article cited above, and it was independently 
10253  used by J. W. Carr III at the University of Michigan 
10264  in 1953. Several years later, the MANIAC III 
10272  computer was designed to include both kinds of 
10280  arithmetic in its hardware; see R. L. Ashenhurst 
10288  and N. Metropolis, |εJACM |≡6 (1959), 415<428; 
10295  IEEE Transactions on Electronic Computers |π|≡E|≡C<|≡1|≡2 
10301  (1963), 896<901; R. L. Ashenhurst, |εProc. Spring 
10308  Joint Computer Conf. |≡2|≡1 (1962), 195<202. 
10314  |πSee also H. L. Gray and C. Harrison, Jr., |εProc. 
10324  Eastern Joint Computer Conf. |≡1|≡6 (1959), 244<248, 
10331  |πand W. G. Wadey, |εJACM |≡7 (1960), 129<139, 
10339  |πfor further early discussions of unnormalized 
10345  arithmetic.|'!|4|4|4|4For early developments 
10349  in interval arithmetic, and some modi_cations, 
10355  see A. Gibb, |εCACM |≡4 (1961), 319<320; |πB. 
10363  A. Chartres, |εJACM |≡1|≡3 (1966), 386<403; |πand 
10370  the book |εInterval Analysis |πby Ramon E. Moore 
10378  (Prentice-Hall, 1966).|'!|4|4|4|4More recent 
10382  work on ⊗oating-point accuracy is summarized 
10388  in two important papers which can be especially 
10396  recommended for further study: W. M. Kahan, |εProc. 
10404  IFIP Congress (1971), |≡2|≡, 1214<1239; |πR. 
10410  P. Brent, |εIEEE Trans. |π|≡C|≡-|≡2|≡2 (1973), 
10416  601<607. Both papers include useful theory and 
10423  demonstrate that it pays o= in practice.|'{A24}*?*?|∨E|∨X|∨E|∨
10430  R|∨C|∨I|∨S|∨E|'{A12}{H9L11M29}Normalized ⊗oating-point 
10433  arithmetic is assumed unless the contrary is 
10440  speci_ed.)|'{A3}|9|1|≡1|≡.|9|4[|εM|*/|↔O|↔l|\] 
10442  |πProve that identity (7) is a consequence of 
10450  (2) through (6).|'{A3}|9|1|≡2|≡.|9|4[|ε|*/M|↔P|↔c|\|π] 
10454  Use identities (2) through (8) to prove that 
10462  |ε(u|4|↔V|4x)|4|↔V|4(v|4|↔V|4y)|4|¬R|4u|4|↔V|4v 
10463  |πwhenever |εx|4|¬R|40 |πand |εy|4|¬R|40.|'{A3}|9|1|≡3|≡.|9|
10467  4[M|*/|↔P|↔c|\|π] Find eitht-digit ⊗oating-decimal 
10471  numbers |εu, v, |πand |εw |πsuch that|'{A9}|εu|4|↔N|4(v|4|↔N
10478  |4w)|4|=|↔6α=↓|4(u|4|↔N|4v)|4|↔N|4w,|;{A9}|πand 
10480  no exponent over⊗ow or under⊗ow occurs during 
10487  any of these computations.|'{A3}|9|1|≡4|≡.|9|4[|ε|*/|↔O|↔c|\|
10491  π] It is possible to have ⊗oating-point numbers 
10499  |εu, v, |πand |εw |πfor which exponent over⊗ow 
10507  occurs during the calculation of |εu|4|↔N|4(v|4|↔N|4w) 
10513  |πbut not during the calculation of (|εu|4|↔N|4v)|4|↔N|4w?|'
10520  {A3}|9|1|≡6|≡.|9|4|ε[M|*/|↔P|↔P|\] |πAre either 
10523  of the following two identities valid for all 
10531  ⊗oating-point numbers |εu? |π(a) |ε0|4|↔B|4(0|4|↔B|4u)|4α=↓|
10535  4u. |π(b) 1|4|↔n|4(1|4|↔n|4|εu)|4α=↓|4u.|'{A3}|9|1|≡7|≡.|9|4
10538  [|ε|*/M|↔P|↔O|\|π] Let |εu|=|g|*2`|g|≤ |πstand 
10542  for |εu|4|↔N|4u. |πFind ⊗oating-binary numbers 
10547  |εu |πand |εv |πsuch that |ε2(u|=|g|*2`|g|≤|4|↔V|4v|=|g|*2`|g2
10552  )|4|¬W|4(u|4|↔V|4v)|=|g|*2`|g2.|'{A3}|9|1|≡8|≡.|9|4[|*/|↔P|↔c|
10553  \|π] Let |ε|≤e|4α=↓|40.0001; |πwhich of the relations 
10560  |εu|4|"2|4v|9(|≤e), u|4|"2|4v|9(|≤e), u|4|¬V|4|9(|≤e) 
10563  |πhold for the following pairs of base 10, excess 
10572  0, eight-digit ⊗oating-point numbers?|'{A3}!!|1|1|1a)|9|εu|4
10576  α=↓|4(|9|11,|4|→α+↓.31415927),|9v|4α=↓|4(|9|11,|4|→α+↓.31416
10576  000);|'!!|1|1b)|9u|4α=↓|4(|9|10,|4|→α+↓.99997000),|9v|4α=↓|4
10577  (|9|11,|4|→α+↓.10000039);|'!!|4c)|9u|4α=↓|4(24,|4|→α+↓.60225
10578  200),|4|4|1v|4α=↓|4(27,|4|→α+↓.00060225);|'!!|1|1d)|9u|4α=↓|
10579  4(24,|4|→α+↓.60225200),|9v|4α=↓|4(31,|4|→α+↓.00000006);|'
10580  !!|4e)|9u|4α=↓|4(24,|4|→α+↓.60225200),|9v|4α=↓|4(32,|4|→α+↓.
10580  00000000).|'{A3}|9|1|≡9|≡.|9|4[M|*/|↔P|↔P|\] |πProve 
10583  (33) and explain why the conclusion cannot be 
10591  strngthened to |εu|4|¬V|4w(|≤e|β1|4α+↓|4|≤e|β2).|'
folio 313 galley 13
10594  |H*?{U0}{H9L11M29}58320#Computer Programming!(Knuth/A.-W.)!ch
10595  . 4!f. 313!g. 13a|'{A12}|≡1|≡0|≡.|9|4[|εm|*/|↔P|↔C|\] 
10600  |π(W. Kahan.) A certain computer performs ⊗oating-point 
10607  arithmetic without proper rounding, and, in fact, 
10614  its ⊗oating-point multiplication routine ignores 
10619  all but the _rst |εp |πmost signi_cant digits 
10627  of the |ε2p-|πdigit product |εf|βuf|βv. (|πThus 
10633  when |εf|βuf|βv|4|¬W|41/b, |πthe least-signi_cant 
10637  digit of |εu|4|↔N|4v |πalways comes out to be 
10645  zero, due to subsequent normalization.) Show 
10651  that this causes the monotonicity of multiplication 
10658  to fail; i.e., there are positive normalized 
10665  ⊗oating-point numbers |εu, v, w |πsuch that |εu|4|¬W|4v 
10673  |πbut |εu|4|↔N|4w|4|¬Q|4v|4|↔N|4w.|'{A3}|≡1|≡1|≡.|9|4[M|*/|↔P
10675  |↔c|\] |πProve Lemma T.|'{A3}|≡1|≡2|≡.|9|4[|εM|*/|↔P|↔M|\] 
10680  |πCarry out the proof of Theorem B and (46) when 
10690  |¬G|εe|βu|4α_↓|4e|βv|¬G|4|¬R|4p.|'{A3}|π|≡1|≡3|≡.|9|4[|εM|*/|
10691  ↔P|↔C|\] |πSome programming languages (and even 
10697  some computers) make use of ⊗oating-point arithmetic 
10704  only, with no provision for exact calculations 
10711  with integers. If operations on integers are 
10718  desired, we can, of course, represent an integer 
10726  as a ⊗oating-point number; and when the ⊗oating-point 
10734  operations satisfy the basic de_nitions in (9), 
10741  we know that all ⊗oating-point operations will 
10748  be exact,{A3}{H9L11M29}|≡1|≡7|≡.|9|4[|ε|*/|↔P|↔l|\] 
10750  |πWrite a |¬m|¬i|¬x subroutine, |¬f|¬c|¬m|¬p, 
10755  which compares the ⊗oating-point number |εu |πin 
10762  location |¬a|¬c|¬c with the ⊗oating-point number 
10768  |εv |πin register A, and which sets the comparison 
10777  indicator to |¬l|¬e|¬s|¬s|¬, |¬e|¬q|¬u|¬a|¬l|¬, 
10781  or |¬g|¬r|¬e|¬a|¬t|¬e|¬r|¬, according as |εu|4|"2|4v, 
10786  u|4|¬Z|4v, |πor |εu|4|"2|4v (|≤e); |πhere |ε|≤e 
10792  |πis stored in location |¬e|¬p|¬s|¬i|¬l|¬o|¬n 
10797  as a nonnegative _xed-point quantity with the 
10804  decimal point assumed at the left of the word. 
10813  Assume normalized inputs.|'{A3}|≡1|≡8|≡.|9|4[|ε|*/M|↔M|↔c|\|π
10816  ] In unnormalized arithmetic is there a suitable 
10824  number |ε|≤e |πsuch that|'{A9}|εu|4|↔N|4(v|4|↔V|4w)|4|¬V|4(u
10828  |4|↔N|4v)|4|↔V|4(u|4|↔N|4w)!!(|≤e)?|;{A9}|≡1|≡9|≡.|9|4[|*/M|↔
10829  L|↔c|\|π] (W. M. Kahan.) Consider the following 
10836  procedure for ⊗oating-point summation of |εx|β1,|4.|4.|4.|4,
10841  |4x|βn:|'{A9}s|β0|4α=↓|4c|β0|4α=↓|40;!y|βk|4α=↓|4x|βk|4|↔B|4
10842  c|βk|βα_↓|β1, s|βk|4α=↓|4s|βk|βα_↓|β1|4|↔V|4y|βk, 
10844  c|βk|4α=↓|4(s|βk|4|↔B|4s|βk|βα_↓|β1)|4|↔B|4y|βk,!!|πfor!!|ε1
10844  |4|¬E|4k|4|¬E|4n.|;{A9}|πLet the relative errors 
10849  in these operations be de_ned by the equations|'
10857  |εy|βk|4α=↓|4(x|βk|4α_↓|4c|βk|βα_↓|β1)(1|4α+↓|4|≤h|βk),!s|βk
10857  |4α=↓|4(s|βk|βα_↓|β1|4α+↓|4y|βk)(1|4α+↓|4|≤s|βk),!c|βk|4α=↓|
10857  4((s|βk|4α_↓|4s|βk|βα_↓|β1)(1|4α+↓|4|≤g|βk)|4α_↓|4y|βk)(1|4α
10857  +↓|4|≤d|βk),|;{H9L11M29}{A9}|πwhere |ε|¬G|≤g|βk|¬G, 
10860  |¬G|≤h|βk|¬G, |¬G|≤s|βk|¬G|4|¬E|4|≤e. |πProve 
10863  that |εs|βn|4α=↓|4|¬K|ur|)1|¬Ek|¬En|) (1|4α+↓|4|≤u|βk)x|βk, 
10866  |πwhere |ε|¬G|≤u|βk|¬G|4|¬E|42|≤e|4α+↓|4O(n|≤e|g2). 
10868  |π[Theorem C says that if |εb|4α=↓|42 |πand |ε|¬Gs|βk|βα_↓|β
10875  1|¬G|4|¬R|4|¬Gy|βk|¬G |πwe have |εs|βk|βα_↓|β1|4α+↓|4y|βk|4α
10878  =↓|4s|βk|4α_↓|4c|βk |πexactly. But in this exercise 
10884  we want to obtain an estimate which is valid 
10893  |εeven when ⊗oating-point operations are not 
10899  carefully rounded, |πassuming only that each 
10905  operation has bounded relative error.]|'{A3}|≡2|≡0|≡.|9|4[|ε
10910  |*/|↔P|↔C|\|π] (S. Linnainmaa.) Find all |εu, 
10916  v |πfor which |¬G|εu|¬G|4|¬R|4|¬Gv|¬G |πand (47) 
10922  fails.|'{A3}|≡2|≡1|≡.|9|4[|ε|*/M|↔L|↔C|\|π] (T. 
10925  J. Dekker.) Theorem C shows how to do exact addition 
10935  of ⊗oating-binary numbers. Explain how to do 
10942  |εexact multiplication: |πExpress the product 
10947  |εuv |πin the form |εw|4α+↓|4w|¬S, |πwhere |εw 
10954  |πand |εw|¬S |πare computed from two given ⊗oating-binary 
10962  numbers |εu |πand |εv, |πusing only the operations 
10970  |↔V, |↔B, and |↔N.|'{A3}|≡2|≡2|≡.|9|4[|ε|*/M|↔L|↔c|\|π] 
10975  Can drift occur in ⊗oating-point multiplication/division? 
10981  Consider the sequence |εx|β0|4α=↓|4u, x|β2|βn|βα+↓|β1|4α=↓|4
10985  x|β2|βn|4|↔N|4v, x|β2|βn|βα+↓|β2|4α=↓|4x|β2|βn|βα+↓|β1|4|↔n|
10986  4v, |πgiven |εu |πand |εv; |πwhat is the largest 
10995  |εk |πsuch that |εx|βk|4|=|↔6α=↓|4x|βk|βα+↓|β2 
10999  |πis possible?|'{A3}|≡2|≡3|≡.|9|4[|ε|*/M|↔P|↔o|\|π] 
11002  Prove or disprove: |εu|4|↔n|4(u|4|πmod|41)|4α=↓|4|"lu|"L 
11006  |πfor all ⊗oating-point |εu.|'{A3}|≡2|≡4|≡.|9|4[|ε|*/M|↔P|↔p|
11010  \|π] (W. M. Kahan.) De_ne appropriate arithmetic 
11017  operations on intervals [|εa, b], |πwhere |εa 
11024  |πand |εb |πare either nonzero ⊗oating-point 
11030  numbers or the special symbols |→α+↓0, |→α_↓0, 
11037  |→α+↓|¬X, |→α_↓|¬X; furthermore |εa|4|¬E|4b, 
11041  |πand |εa|4α=↓|4b |πimplies that |εa |πis _nite 
11048  and nonzero. This interval stands for all ⊗oating-point 
11056  |εx |πsuch that |εa|4|¬E|4x|4|¬E|4b, |πwhere 
11061  we regard |→α_↓|¬X|4|¬W|4|→α_↓|εu|4|¬W|4|→α_↓0|4|¬W|40|4|→α+
11063  ↓0|4|¬W|4|→α+↓u|4|¬W|4|→α+↓|¬X |πfor all positive 
11067  |εu. (|πThus, [1,|42] means |ε1|4|¬E|4x|4|¬E|42, 
11072  [|→α+↓0,|41] |πmeans |ε0|4|¬W|4x|4|¬E|41, [|→α_↓0,|41] 
11076  |πmeans |ε0|4|¬E|4x|4|¬E|41, [|→α_↓0,|4|→α+↓0] 
11079  |πdenotes the single value 0, and [|→α_↓|¬X, 
11086  |→α+↓|¬X] stands for everything.)|'{A3}|≡2|≡5|≡.|9|4[|ε|*/|↔O
11090  |↔C|\|π] When people speak about inaccuracy in 
11097  ⊗oating-point arithmetic they often ascribe errors 
11103  to ``cancellation'' which occurs during the subtraction 
11110  of nearly equal quantities. But when |εu |πand 
11118  |εv |πare approximately equal, the di=erence 
11124  |εu|4|↔B|4v |πis obtained exactly, with no error*3 
11131  What do these people really mean?|'{A3}|≡2|≡6|≡.|9|4[|ε|*/M|↔
11137  P|↔C|\|π] (W. M. Kahan.) Let |εf|1(x)|4α=↓|41|4α+↓|4x|4α+↓|4
11142  x|g2|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4x|g1|g0|g6|4α=↓|4(1|4α_↓|4x|g
11142  1|g0|g7)/(1|4α_↓|4x) |πfor |εx|4|¬W|41, |πand 
11146  let |εg(y)|4α=↓|4f|1{H12}({H10}(|f1|d32|)|4α_↓|4y|g2)(3|4α+↓
11147  |43.45y|g2){H12}) |π{H10}for 0|4|¬W|4|εy|4|¬W|41. 
11150  |πEvaluate |εg(y) |πon one or more ``pocket calculators,'' 
11158  for |εy|4α=↓|410|gα_↓|g3, 10|gα_↓|g4, 10|gα_↓|g5, 
11162  10|gα_↓|g6, |πand explain all inaccuracies in 
11168  the results you obtain. [Since present-day calculators 
11175  do not round correctly, the results are often 
11183  surprising. Note that |εg(|≤e)|4α=↓|4107|4α_↓|410491.35|≤e|g
11186  2|4α+↓|4O(|≤e|g4).]|'|H*?*?{U0}{H10L12M29}58320#Computer 
folio 317 galley 14
11188  Programming!(Knuth/A.-W.)!ch. 4!f. 317!g. 14a|'
11192  |≤⊂|π|∨4|∨.|∨2|∨.|∨3|∨.|9|∨D|∨o|∨u|∨b|∨l|∨e|∨-|∨P|∨r|∨e|∨c|∨
11192  i|∨s|∨i|∨o|∨n|9|∨C|∨a|∨l|∨c|∨u|∨l|∨a|∨t|∨i|∨o|∨n|∨s|'
11193  {A6}Up to now we have considered ``single-precision'' 
11200  ⊗oating-point arithmetic, which essentially means 
11205  that the ⊗oating-point values we have dealt with 
11213  can be stored in a single machine word. When 
11222  single-precision ⊗oating-point arithmetic does 
11226  not yield su∃cient accuracy for a given application, 
11234  the precision can be increased by suitable programming 
11242  techniques which use two or more words of memory 
11251  to represent each number. Although we shall discuss 
11259  the general question of high-precision calculations 
11265  in Section 4.3, it is appropriate to give a separate 
11275  discussion of double-precision here; special 
11280  techniques apply to double precision which are 
11287  comparatively inappropriate for higher precisions, 
11292  and double precision is a reasonably important 
11299  topic in its own right, since it is the _rst 
11309  step beyond single precision and is applicable 
11316  to many problems which do not require extremely 
11324  high precision.|'!|4|4|4|4Double-precision calculations 
11328  are almost always required for ⊗oating-point, 
11334  rather than _xed-point arithmetic, except perhaps 
11340  in statistical work where _xed-point double-precision 
11346  arithmetic is commonly used to calculate sums 
11353  of squares and cross products; since _xed-point 
11360  versions of double-precision arithmetic are simpler 
11366  than ⊗oating-point versions, we will con_ne our 
11373  discussion here to the latter.|'!|4|4|4|4Double 
11379  precision is quite frequently desired not only 
11386  to extend the precision of the fraction parts 
11394  of ⊗oating-point numbers, but also to increase 
11401  the range of the exponent part. Thus we shall 
11410  deal in this section with the following two-word 
11418  format for double-precision ⊗oating-point numbers 
11423  in the |¬m|¬i|¬x computer:|'{A9}|¬N|ε!!e!!e!!f!!f!!f!!!!f!!f
11427  !!f!!f!!f!|9.|J!(1)|;{A12}|πHere two bytes are 
11432  used for the exponent and eight bytes are used 
11441  for the fraction. The exponent is ``excess |εb|g2/2,'' 
11449  |πwhere |εb |πis the byte size. The sign will 
11458  appear in the most signi_cant word; it is convenient 
11467  to ignore the sign of the other word completely.|'
11476  !|4|4|4Our discussion of double-precision arithmetic 
11481  will be rather machine-oriented, because it is 
11488  only by studying the problems involved in coding 
11496  these routines that a person can properly appreciate 
11504  the subject. A careful study of the |¬m|¬i|¬x 
11512  programs below is therefore essential to the 
11519  understanding of this section.|'!|4|4|4|4In this 
11525  section we shall depart from the idealistic goals 
11533  of accuracy stated in the previous two sections; 
11541  our double-precision routines will |εnot |πround 
11547  their results, and a little bit of error will 
11556  sometimes be allowed to creep in*3 Users dare 
11564  not trust these routines too much. There was 
11572  ample reason to squeeze out every possible drop 
11580  of accuracy in the single-precision case, but 
11587  now we face a di=erent situation: (a) The extra 
11596  programming required to ensure true double-precision 
11602  rounding in all cases is considerable; fully 
11609  accurate routines would take, say, twice as much 
11617  space and half again as much more time. It was 
11627  comparatively eacy to make our single-precision 
11633  routines perfects, but double precision brings 
11639  us face to face with our machine's limitations. 
11647  [A similar situation occurs with respect to other 
11655  ⊗oating-point subroutines; we can't expect the 
11661  routine to compute round(cos|4|εx) |πexactly 
11666  for all |εx, |πsince that turns out to be virtually 
11676  impossible. Instead, the cosine routine should 
11682  provide the best relative error it can achieve 
11690  with reasonable speed, for all reasonable values 
11697  of |εx; |πand the designer of the routine should 
11706  try to make the computed function satisfy simple 
11714  mathematical laws whenever possible (e.g., cos 
11720  (|→α_↓|εx)|4α=↓|4|πcos|4|εx, |¬G|πcos|4|εx|¬G|4|¬E|41, 
11722  |πcos x|4|¬R|4cos|4|εy |πfor |ε0|4|¬E|4x|4|¬E|4y|4|¬W|4|≤p).
11725  ] |π(b) Single-precision arithmetic is a ``staple 
11732  food'' that everybody who wants to employ ⊗oating-point 
11740  arithmetic must use, but double precision is 
11747  usually for situations where such clean results 
11754  aren't as important. The di=erence between seven- 
11761  and eight-place accuracy can be noticeable, but 
11768  we rarely care about the di=erence between 15- 
11776  and 16-place accuracy. Double precision is most 
11783  often used for intermediate steps during the 
11790  calculation of single-precision results; its 
11795  full potential isn't needed. (c) It will be instructive 
11804  for us to analyze these procedures in order to 
11813  see how inaccurate they can be, since they typify 
11822  the types of shor cuts generally taken in bad 
11831  single-precision routines (see exercises 7 and 
11837  8).|'!|4|4|4|4Let us now consider addition and 
11844  subtraction operations from this standpoint. 
11849  Double-precision addition di=ers from the single-precision 
11855  case in the following respects: The addition 
11862  is performed by separately adding together the 
11869  least-signi_cant halves and the most-signi_cant 
11874  halves, taking care of carries appropriately. 
11880  But since we are doing signed-magnitude arithmetic, 
11887  it is possible to add the least-signi_cant halves 
11895  and to get the wrong sign (namely, when the signs 
11905  of the operands are opposite, and the least-signi_cant 
11913  half of the smaller operand is bigger than the 
11922  least-signi_cant half of the larger operand); 
11928  so it is helpful to know what sign to expect. 
11938  Therefore in step A2 (cf. Algorithm 4.2.1A), 
11945  we not only assume that |εe|βu|4|¬R|4e|βv, |πwe 
11952  also assume that |¬G|εu|¬G|4|¬R|4|¬Gv|¬G. |πThis 
11957  means we can be sure that the _nal sign will 
11967  be the sign of |εu. |πIn other respects, the 
11976  double-precision addition is very much like the 
11983  single-precision routine, only everything is 
11988  done twice.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡A (|εDouble-precis
11992  ion addition). |πThe subroutine |¬d|¬f|¬a|¬d|¬d 
11997  adds a double-precision ⊗oating-point number 
12002  |εv, |πhaving the form (1), to a double-precision 
12010  ⊗oating-point number |εu, |πassuming that |εv 
12016  |πis initially in rAX (i.e., registers A and 
12024  X), and that |εu |πis initially stored in locations 
12033  |¬a|¬c|¬c and |¬a|¬c|¬c|¬x. The answer appears 
12039  both in rAX and in (|¬a|¬c|¬c, |¬a|¬c|¬c|¬x). 
12046  The subroutine |¬d|¬f|¬s|¬u|¬b subtracts |εv 
12051  |πfrom |εu |πunder the same conventions. Both 
12058  input operands are assumed to be normalized, 
12065  and the answer is normalized. The last portion 
12073  of this program is a double-precision normalization 
12080  procedure which is used by other subroutines 
12087  of this section. Exercise 5 shows how to improve 
12096  this program signi_cantly.|'{A12}{H9L11M29}|∂!!|∂!!!!!|∂!!!!
12099  |∂!!!!!!!|∂!!!!!!!!!!!!!!!!!!!!|4|4|4|∂|E|;|>
12101  |ε|*/|↔c|↔O|'|\|¬a|¬b|¬s|'|¬e|¬q|¬u|'|¬1|¬.|¬5|'
12105  |πField|4de_nition|4for|4absolute|4value|'>|>
12108  |ε|*/|↔c|↔P|'|π|\|¬s|¬i|¬g|¬n|'|¬e|¬q|¬u|'|¬0|¬.|¬0|'
12112  Field|4de_nition|4for|4sign|'>|>|ε|*/|↔c|↔L|'|π|\|¬e|¬x|¬p|¬d
12116  |'|¬e|¬q|¬u|'|¬1|¬.|¬2|'Double-precision exponent 
12121  _eld|'>|>|ε|*/|↔c|↔M|'|π|\|¬d|¬f|¬s|¬u|¬b|'|¬s|¬t|¬a|'
12127  |¬t|¬e|¬m|¬p|'Double-precision subtraction:|'
12130  >|>|ε|*/|↔c|↔C|'|π|\|;|¬l|¬d|¬a|¬n|'|¬t|¬e|¬m|¬p|'
12136  Change sign of |εv.|'>|>|ε|*/|↔c|↔o|'|π|\|¬d|¬f|¬a|¬d|¬d|'
12144  |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'Double-precision 
12147  addition:|'>|>|ε|*/|↔c|↔p|'|π|\|;|¬c|¬m|¬p|¬a|'
12153  |¬a|¬c|¬c|≤#|¬a|¬b|¬s|≤&|'Compare |ε|¬Gu|¬G |πwith 
12157  |ε|¬Gv|¬G.|'>|>|*/|↔c|↔l|'|\|π|;|¬j|¬g|'|¬1|¬f|'
12164  >|>|ε|*/|↔c|↔m|'|\|π|;|¬j|¬l|'|¬2|¬f|'>|>|ε|*/|↔O|↔c|'
12173  |\|π|;|¬c|¬m|¬p|¬x|'|¬a|¬c|¬c|¬x|≤#|¬a|¬b|¬s|≤&|'
12176  >|>|ε|*/|↔O|↔O|'|\|π|;|¬j|¬l|¬e|'|¬2|¬f|'>|>|ε|*/|↔O|↔P|'
12185  |\|π|¬1|¬h|'|¬s|¬t|¬a|'|¬a|¬r|¬g|'If |¬G|εu|¬G|4|¬W|4|¬Gv|¬G
12189  , |πinterchange|4|εu|4|"m|4v.|'>|>|ε|*/|↔O|↔L|'
12194  |π|\|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'>|>|ε|*/|↔O|↔M|'
12200  |π|\|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔O|↔C|'|\|π|;
12207  |¬l|¬d|¬x|'|¬a|¬c|¬c|¬x|'>|>|ε|*/|↔O|↔o|'|\|π|;
12213  |¬e|¬n|¬t|¬1|'|¬a|¬c|¬c|'(|¬a|¬c|¬c and|4|¬a|¬c|¬c|¬x|4are|4
12216  in|4consecutive|'>|>|ε|*/|↔O|↔p|'|\|π|;|¬m|¬o|¬v|¬e|'
12222  |¬a|¬r|¬g|≤#|¬2|≤&|'!!locations.)|'>|>|ε!|↔O|↔p|'
12227  |\|π|¬2|¬h|'|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'Now|4|¬a|¬c|¬c|4has|4th
12230  e|4sign|4of|4the|4answer.|'>|>|ε|*/|↔O|↔m|'|\|π|;
12235  |¬l|¬d|¬1|¬n|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|¬d|≤&|'
12237  rI1|4|¬L|4|→α_↓|εe|βv.|'>|>|ε|*/|↔P|↔c|'|\|π|;
12242  |¬l|¬d|¬2|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|¬d|≤&|'rI2|4|¬L|4|εe|βu.|'
12245  >|>|ε|*/|↔P|↔O|'|\|π|;|¬i|¬n|¬c|¬1|'|¬0|¬,|¬2|'
12251  rI1|4|¬L|4e|βu|4α_↓|4e|βv.|'>|>|ε|*/|↔P|↔P|'|\|π|;
12256  |¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'>|>|ε|*/|↔P|↔L|'
12262  |\|π|;|¬s|¬r|¬a|¬x|'|¬1|¬,|¬1|'Scale|4right|'
12266  >|>|ε|*/|↔P|↔M|'|\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|'0|4|εv|β1|4v|β2|4
12272  v|β3|4v|β4|'>|>|ε|*/|↔P|↔C|'|\|π|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'
12279  |εv|β5|4v|β6|4v|β7|4v|β8|4v|β9|'>|>|ε|*/|↔P|↔o|'
12283  |\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|¬x|≤#|¬s|¬i|¬g|¬n|≤&|'
12286  Store|4true|4sign|4in|4both|4halves.|'>|>|ε|*/|↔P|↔p|'
12290  |\|π|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔P|↔l|'|\|π|;
12297  |¬l|¬d|¬x|'|¬a|¬c|¬c|¬x|'>|>|ε|*/|↔P|↔m|'|\|π|;
12303  |¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'>*?*?*?|>|ε|*/|↔L|↔c|'
12309  |\|π|;|¬s|¬t|¬a|'|¬a|¬c|¬c|'|εu|β1|4u|β2|4u|β3|4u|β4|4u|β5|'
12313  >|>|ε|*/|↔L|↔O|'|;|\|¬s|¬l|¬a|¬x|'|¬4|'>|>|ε|*/|↔L|↔P|'
12322  |\|;|¬e|¬n|¬t|¬x|'|¬1|'>|>|ε|*/|↔L|↔L|'|\|;|¬s|¬t|¬x|'
12330  |¬e|¬x|¬p|¬o|'|¬e|¬x|¬p|¬o|4|¬L|41|4|π(see below).|'
12333  >|>|ε|*/|↔L|↔M|'|\|;|¬s|¬r|¬c|'|¬1|'1|4u|β5|4u|β6|4u|β7|4u|β8
12339  |'>|>|ε|*/|↔L|↔C|'|\|;|¬s|¬t|¬a|'|¬1|¬f|≤#|¬s|¬i|¬g|¬n|≤&|'
12346  |πA|4trick,|4see|4comments|4in|4text.|'>|>|ε|*/|↔L|↔o|'
12350  |\|;|¬a|¬d|¬d|'|¬a|¬r|¬g|¬x|≤#|¬0|¬.|¬4|≤&|'|πAdd|4|ε0|4v|β5
12353  |4v|β6|4v|β7|4v|β8.|'>|>|ε|*/|↔L|↔p|'|\|;|¬s|¬r|¬a|¬x|'
12359  |¬4|'>|>|ε|*/|↔L|↔l|'|\|¬1|¬h|'|¬d|¬e|¬c|¬a|'|¬1|'
12366  |πRecover|4from|4inserted|41|4(Sign|4varies)|'
12367  >|>|ε|*/|↔L|↔m|'|\|;|¬a|¬d|¬d|'|¬a|¬c|¬c|≤#|¬0|¬.|¬4|≤&|'
12373  |πAdd|4most|4signi_cant|4halves.|'>|>|ε|*/|↔M|↔c|'
12377  |\|π|;|¬a|¬d|¬d|'|¬a|¬r|¬g|'(Over⊗ow|4cannot|4occur)|'
12381  >|>|ε|*/|↔M|↔O|'|\|π|¬d|¬n|¬o|¬r|¬m|'|¬j|¬a|¬n|¬z|'
12386  |¬1|¬f|'Normalization|4routine:|'>|>|ε|*/|↔M|↔P|'
12391  |\|π|;|¬j|¬x|¬n|¬z|'|¬1|¬f|'|εf|βw|4|πin|4rAX,|4|εe|βw|4α=↓|
12394  4|¬e|¬x|¬p|¬o|4α+↓|4rI2.|'>|>|ε|*/|↔M|↔L|'|\|π|¬d|¬z|¬e|¬r|¬o
12398  |'|¬s|¬t|¬a|'|¬a|¬c|¬c|'If|4|εf|βw|4α=↓|40,|4|πset|4|εe|βw|4
12401  |¬L|40.|'>|>|ε|*/|↔M|↔M|'|\|π|;|¬j|¬m|¬p|'|¬9|¬f|'
12408  >|>|ε|*/|↔M|↔C|'|\|π|¬2|¬h|'|¬s|¬l|¬a|¬x|'|¬1|'
12414  Normalize|4to|4left.|'>|>|ε|*/|↔M|↔o|'|\|π|;|¬d|¬e|¬c|¬2|'
12420  |¬1|'>|>|ε|*/|↔M|↔p|'|\|π|¬1|¬h|'|¬c|¬m|¬p|¬a|'
12426  |→|≤$|¬0|≤$|≤#|¬1|¬.|¬1|≤&|'Is|4the|4leading|4byte|4zero?|'
12428  >|>|ε|*/|↔M|↔l|'|\|π|;|¬j|¬e|'|¬2|¬b|'>|>|ε|*/|↔M|↔m|'
12437  |\|π|;|¬j|¬e|'|¬2|¬b>|>|ε|*/|↔M|↔m|'|\|π|;|¬s|¬r|¬a|¬x|'
12444  |¬2|'(Rounding|4omitted)|'>|>|ε|*/|↔C|↔c|'|\|π|;
12450  |¬s|¬t|¬a|'|¬a|¬c|¬c|'>|>|ε|*/|↔C|↔O|'|\|π|;|¬l|¬d|¬a|'
12457  |¬e|¬x|¬p|¬o|'Compute|4_nal|4exponent.|'>|>|ε|*/|↔C|↔P|'
12462  |\|π|;|¬i|¬n|¬c|¬a|'|¬0|¬,|¬2|'>|>|ε|*/|↔C|↔L|'
12468  |\|π|;|¬j|¬a|¬n|'|¬e|¬x|¬p|¬u|¬n|¬d|'Is|4it|4negative?|'
12472  >|>|ε|*/|↔C|↔M|'|\|π|;|¬s|¬t|¬a|'|¬a|¬c|¬c|≤#|¬e|¬x|¬p|¬d|≤&|
12477  '>|>|ε|*/|↔C|↔C|'|\|ε|;|¬c|¬m|¬p|¬a|'|→|≤$|¬1|≤#|¬3|¬.|¬3|≤&|
12483  →|≤$|'|πIs|4it|4more|4than|4two|4bytes?|'>|>|ε|*/|↔C|↔o|'
12488  |\|π|;|¬j|¬l|'|¬8|¬f|'>|>|ε|*/|↔C|↔p|'|\|π|¬e|¬x|¬p|¬o|¬v|¬d|
12494  '|¬h|¬l|¬t|'|¬2|¬0|'>|>|ε|*/|↔C|↔l|'|\|π|¬e|¬x|¬p|¬u|¬n|¬d|'
12501  |¬h|¬l|¬t|'|¬1|¬0|'>|>|ε|*/|↔C|↔m|'|\|π|¬8|¬h|'
12507  |¬l|¬d|¬a|'|¬a|¬c|¬c|'Bring|4answer|4into|4rA.|'
12510  >|>|ε|*/|↔o|↔c|'|\|π|¬9|¬h|'|¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'
12516  >|>|ε|*/|↔o|↔O|'|\|π|¬e|¬x|¬i|¬t|¬d|¬f|'|¬j|¬m|¬pα/↓|'
12521  Exit|4from|4subroutine.|'>|>|ε|*/|↔o|↔P|'|\|π|¬a|¬r|¬g|'
12526  |¬c|¬o|¬n|'|¬0|'>|>|ε|*/|↔o|↔L|'|\|π|¬a|¬r|¬g|¬x|'
12532  |¬c|¬o|¬n|'|¬0|'>|>|ε|*/|↔o|↔M|'|\|π|¬a|¬c|¬c|'
12538  |¬c|¬o|¬n|'|¬0|'Floating-point|4accumulator|'
12541  >|>|ε|*/|↔o|↔C|'|\|π|¬a|¬c|¬x|'|¬c|¬o|¬n|'|¬0|'
12547  >|>|ε|*/|↔o|↔o|'|\|π|¬e|¬x|¬p|¬o|'|¬c|¬o|¬n|'|¬0Part|4of|4``r
12552  aw|4exponent''|'>|H*?*?*?*?{U0}{H10L12M29}58320#Computer 
folio 322 galley 15-16 WARNING: End of 15 and beginning of 16 were unreadable.
12555  Programming!(Knuth/A.-W.)!ch. 4!f. 322!g. 15a|'
12559  {A12}!|4|4|4|4When the least-signi_cant halves 
12563  are added together in this program, an extra 
12571  digit ``1'' is inserted at the left of the word 
12581  which is known to have the correct sign. After 
12590  the addition, this byte can be 0, 1, or 2, depending 
12601  on the circumstances, and all three cases are 
12609  handled simultaneously in this way. (Compare 
12615  this with the rather cumbersome method of complementation 
12623  which is used in Program 4.2.1A.)|'!|4|4|4|4It 
12630  is worth noting that register A can be zero after 
12640  the instruction on line 40 has been performed; 
12648  and, because of the way |¬m|¬i|¬x de_nes the 
12656  sign of a zero result, the accumulator contains 
12664  the correct sign which is to be attached to the 
12674  result if register X is nonzero. If lines 39 
12683  and 40 were interchanged, the program would be 
12691  incorrect (although both instructions are ``|¬a|¬d|¬d'')*3|'
12697  !|4|4|4|4Now let us consider double-precision 
12702  multiplication. The product has four components, 
12708  shown schematically in Fig. 4. Since we need 
12716  only the leftmost eight bytes, it is convenient 
12724  to work only with the digits to the left of the 
12735  vertical line in the diagram, and this means 
12743  in particular that we need not even compute the 
12752  product of the two least signi_cant halves.|'
12759  {A6}{H9L11M29}|≡F|≡i|≡g|≡.|4|1|≡4|≡.|9|4Double-precision 
12760  multiplication of witht1byte fraction parts.|;
12765  {A12}{H10L12M29}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡M (|εDouble-precisio
12767  n multiplication).|9|4|πThe input and output 
12772  conventions for this subroutine are the same 
12779  as for Program A.|'{A12}{H9L11M29}|∂!!|∂!!!!|∂!!!!|∂!!!!!!!!
12783  !|∂!!!!!!!!!!!!!!!!!!!|4|4|4|∂|E|;|>|ε|*/|↔c|↔O|'
12786  |\|π|¬b|¬y|¬t|¬e|'|¬e|¬q|¬u|'|¬1|≤#|¬4|¬.|¬4|≤&|'
12789  Byte|4size|'>|>|ε|*/|↔c|↔P|'|\|π|¬q|¬q|'|¬e|¬q|¬u|'
12795  |¬b|¬y|¬t|¬e{J3}|≤∩|¬b|¬y|¬t|¬e/|¬2|'Excess|4of|4double|4pre
12796  cision|4exponent|'>|>|ε|*/|↔c|↔L|'|\|π|¬d|¬f|¬m|¬u|¬l|'
12801  |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'Double-precision|4multiplicat
12803  ion:|'>|>|ε|*/|↔c|↔M|'|\|π|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'
12810  >|>|ε|*/|↔c|↔C|'|\|π|;|¬s|¬l|¬a|¬x|'|¬2|'Remove|4exponent.|'
12817  |>|ε|*/|↔c|↔o|'|\|π|;|¬s|¬t|¬a|'|¬a|¬r|¬g|'|εv|βm|'
12823  >|>|*/|↔c|↔l|'|\|π|;|¬s|¬t|¬x|'|¬a|¬r|¬g|¬x|'|εv|βl|'
12830  >|>|*/|↔c|↔l|'|\|;|¬l|¬d|¬a|'|¬t|¬e|¬m|¬p|≤#|¬e|¬x|¬p|¬d|≤&|'
12836  >|>|*/|↔C|↔M|'|\|;|¬A|¬D|¬D|'|¬A|¬C|¬C|≤#|¬E|¬X|¬P|¬D|≤&|'
12842  >|>|*/|↔O|↔c|'|\|;|¬s|¬t|¬a|'|¬e|¬x|¬p|¬o|'|¬e|¬x|¬p|¬o|4|¬L|
12848  4e|βu|4α+↓|4e|βv.|'>|>|*/|↔O|↔O|'|\|;|¬e|¬n|¬t|¬2|'
12854  |→|≤*/|¬q|¬q|'|πrI2|4|¬L|4|→α_↓|¬q|¬q|'>|>|*/|↔O|↔P|'
12859  |\|;|¬l|¬d|¬a|'|¬a|¬c|¬c|'>|>|*/|↔O|↔L|'|\|;|¬l|¬d|¬x|'
12867  |¬a|¬c|¬c|¬x|'>|>|ε|*/|↔O|↔M|'|\|;|¬s|¬l|¬a|¬x|'
12873  |¬2|'|πRemove|4exponent.|'>|>|ε|*/|↔O|↔C|'|\|;
12879  |¬s|¬t|¬a|'|¬a|¬c|¬c|'u|βm|'>|>|ε|*/|↔O|↔o|'|\|;
12886  |¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'u|βl|'>|>|ε|*/|↔O|↔p|'
12892  |\|;|¬m|¬u|¬l|'|¬a|¬r|¬g|¬x|'u|βm|4α⊗↓|4v|βl|'
12896  >|>|ε|*/|↔O|↔l|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'>
12903  |>|ε|*/|↔O|↔m|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|≤#|¬a|¬b|¬s|≤&|'
12908  >|>|ε|*/|↔P|↔c|'|\|;|¬m|¬u|¬l|'|¬a|¬c|¬c|¬x|≤#|¬a|¬b|¬s|≤&|'
12914  |¬Gv|βm|4α⊗↓|4u|βl|¬G|'>|>|ε|*/|↔P|↔O|'|\|;|¬s|¬r|¬a|'
12920  |¬1|'0|4x|4x|4x|4x|'>|>|ε|*/|↔P|↔P|'|;|\|¬a|¬d|¬d|'
12927  |¬t|¬e|¬m|¬p|≤#|¬1|¬.|¬4|≤&|'(|πOver⊗ow|4cannot|4occur)|'
12929  >|>|ε|*/|↔P|↔L|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|'>
12936  |>|ε|*/|↔P|↔M|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|'>|>|ε|*/|↔P|↔C|'
12944  |\|;|¬m|¬u|¬l|'|¬a|¬c|¬c|'v|βm|4α⊗↓|4u|βm|'>|>
12950  |ε|*/|↔P|↔o|'|\|;|¬s|¬t|¬a|'|¬t|¬e|¬m|¬p|≤#|¬s|¬i|¬g|¬n|≤&|'
12954  |πTrue|4sign|4of|4result.|'>|>|ε|*/|↔P|↔p|'|\|;
12959  |¬s|¬t|¬a|'|¬a|¬c|¬c|'|πNow|4prepare|4to|4add|4all|4the|'
12962  >|>|ε|*/|↔P|↔l|'|\|;|¬s|¬t|¬x|'|¬a|¬c|¬c|¬x|'!!|πpartial|4pro
12968  ducts|4together.|'>|>|ε|*/|↔P|↔m|'|\|;|¬l|¬d|¬a|'
12974  |¬a|¬c|¬c|¬x|≤#|¬0|¬.|¬4|≤&|'0|4x|4x|4x|4x|'>
12977  |>|ε|*/|↔L|↔c|'|\|;|¬a|¬d|¬d|'|¬t|¬e|¬m|¬p|'(|πOver⊗ow|4canno
12982  t|4occur)|'>|>|ε|*/|↔L|↔O|'|\|;|¬s|¬r|¬a|¬x|'|¬4|'
12989  >|>|ε|*/|↔L|↔P|'|\|;|¬a|¬d|¬d|'|¬a|¬c|¬c|'|π(Over⊗ow|4cannot|
12995  4occur)|'>|>|ε|*/|↔L|↔L|'|\|;|¬j|¬m|¬p|'|¬d|¬n|¬o|¬r|¬m|'
13002  |πNormalize|4and|4exit.|'>{A12}{H10L12M29}!|4|4|4|4Note 
13005  the careful treatment of signs in this program, 
13013  and note also the fact that the range of exponents 
13023  makes it impossible to compute the _nal exponent 
13031  using an index register. Program M is perhaps 
13039  a little too slipshod in accuracy, since it throws 
13048  away all the information to the right of the 
13057  vertical line in Fig. 4 and this can make the 
13067  least signi_cant byte as much as 2 in error. 
13076  A little more accuracy can be achieved as discussed 
13085  in exercise 4.|'{A12}!|4|4|4|4Double-precision 
13089  ⊗oating division is the most di∃cult routine, 
13096  or at least the most frightening prospect we 
13104  have encountered so far in this chapter. Actually, 
13112  it is not terribly complicated, once we see how 
13121  to do it; let us write the numbers to be divided 
13132  in the form |ε(u|βm|4α+↓|4|≤eu|βl)/(v|βm|4α+↓|4|≤ev|βl), 
13136  |πwhere |ε|≤e |πis the reciprocal of the word 
13144  size of the computer, and where |εv|βm |πis assumed 
13153  to be normalized. The fraction can now be expanded 
13162  as follows:|'{A9}|h|εu|βm|4α+↓|4|≤eu|βl|4|∂α=↓|4u|βm|4α+↓|4|
13164  ≤eu|βl|4|4|4|4|41|4α_↓|4|≤e|4|4|4v|βm|4|4|4|4α+↓|4|≤e|g2|4|4
13164  |4v|βm|4|4|4|g2|4α_↓|4.|4.|4.|4|4|4.|E|n|;| |(u|βm|4α+↓|4|≤e
13165  u|βl|d2v|βm|4α+↓|4|≤ev|βl|)|4|Lα=↓|4|(u|βm|4α+↓|4|≤eu|βl|d2v
13165  |βm|)|4|↔a|(1|d21|4α+↓|4|≤e(v|βl/v|βm)|)|↔s>{A5}|Lα=↓|4|(u|β
13166  m|4α+↓|4|≤eu|βl|d2v|βm|)|4{H12}|↔a{H10}1|4α_↓|4|≤e|↔a|(v|βl|
13166  d2v|βm|)|↔s|4α+↓|4|≤e|g2|↔a|(v|βl|d2v|βm|)|↔s|g2|4α_↓|4|¬O|4
13166  |¬O|4|¬O{H12}|↔s{H10}.|J!(2)>{A9}|πSince |ε0|4|¬E|4|¬Gv|βl|¬
13168  G|4|¬W|41 |πand |ε1/b|4|¬E|4|¬Gv|βm|¬G|4|¬W|41, 
13171  |πwe have |ε|¬Gv|βl/v|βm|¬G|4|¬W|4b, |πand the 
13176  error from dropping terms in |ε|≤e|g2 |πcan be 
13184  disregarded. Our method therefore is to compute 
13191  |εw|βm|4α+↓|4|≤ew|βl|4α=↓|4(u|βm|4α+↓|4|≤eu|βl)/v|βm, 
13192  |πand then to subtract |ε|≤e |πtimes |εw|βmv|βl/v|βm 
13199  |πfrom the result.|'!|4|4|4|4In the following 
13205  program, lines 27<32 do the lower half of a double-preision 
13215  addition, using another method for forcing the 
13222  appropriate sign as an alternative to the trick 
13230  of Program A.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡D 
13235  (|εDouble-precision division).|4|9|πThis program 
13238  adheres to the same conventions as Programs A 
13246  and M.|'{A12}{H9L11M29}|∂!!|∂!!!!!|∂!!!!|∂!!!!!!!!!!!|∂!!!!!
13248  !!!!!!!!!!!|4|4|4|∂|E|;|>|ε|*/|↔c|↔O|'|\|¬d|¬f|¬d|¬i|¬v|'
13252  |¬s|¬t|¬j|'|¬e|¬x|¬i|¬t|¬d|¬f|'|πDouble-precision|4division:
13254  |'>|>|ε|*/|↔P|↔c|'|\|;|¬l|¬d|¬a|'|¬a|¬r|¬g|¬x|≤#|¬1|¬.|¬4|≤&|
13260  '>*?uble-precision division technique by hand, 
13267  with |ε|≤e|4α=↓|4|f1|d31000|), |πwhen dividing 
13271  180000 by 314159. {H11}({H9}Thus, let (|εu|βm,|4u|βl)|4α=↓|4
13276  (.180,|4.000) |πand |ε(v|βm,|4v|βl)|4α=↓|4(.314,|4.159), 
13279  |πand _nd the quotient using the method suggested 
13287  in the text following (2).{H11})|'{A3}{H9}|9|1|≡2|≡.|9|4[|ε|
13292  */|↔P|↔c|\] |πWould it be a good idea to insert 
13301  the instruction ``|¬e|¬n|¬t|¬x |¬0'' between 
13306  lines 30 and 31 of Program M, in order to keep 
13317  unwanted information left over in register X 
13324  from interfering with the accuracy of the results?|'
13332  {A3}{H9}|9|1|≡3|≡.|9|4[|ε|*/M|↔P|↔c|\] |πExplain 
13334  why over⊗ow cannot occur during Program M.|'{A3}|9|1|≡4|≡.|9
13341  |4[|ε|*/|↔P|↔P|\] |πHow should Program M be changed 
13348  so that extra accuracy is achieved, essentially 
13355  by moving the vertical line in Fig. 4 over to 
13365  the right one position? Specify all changes that 
13373  are required, and determine the di=erence in 
13380  execution time caused by these changes.|'{A3}|9|1|≡5|≡.|9|4[
13386  |ε|*/|↔P|↔M|\] |πHow should Program A be changed 
13393  so that extra accuracy is achieved, essentially 
13400  by working with a nine-byte accumulator instead 
13407  of an eight-byte accumulator to the right of 
13415  the decimal point? Specify all changes that are 
13423  required, and determine the di=erence in execution 
13430  time caused by these changes.|'{A3}{H9}|9|1|≡6|≡.|9|4[|ε|*/|↔
13435  P|↔L|\] |πAssume that the double-precision subroutines 
13441  of this section and the single-precision subroutines 
13448  of Section 4.2.1 are being used in the same main 
13458  program. Write a subroutine which converts a 
13465  single-precision ⊗oating-point number into double-precision 
13470  form (1), and write another subroutine which 
13477  converts a double-precision ⊗oating-point number 
13482  into single-precision form (or which reports 
13488  exponent over⊗ow or under⊗ow if the conversion 
13495  is impossible).|'{A3}{H9}|9|1|≡7|≡.|9|4[|ε|*/M|↔L|↔c|\|π] 
13498  Estimate the accuracy of the double-precision 
13504  subroutines in this section, by _nding bounds 
13511  |ε|≤d|β1, |≤d|β2, |πand |ε|≤d|β3 |πon the relative 
13518  errors|'{A9}|ε|¬G{H11}({H9}(u|4|↔V|4v)|4α_↓|4(u|4α+↓|4v){H11
13519  })|4h9}/(u|4α+↓|4v)|¬G,!!|¬G{H11}({H9}(u|4|↔N|4v)|4α_↓|4(u|4
13519  α⊗↓|4v){H11}){H9}/(u|4α⊗↓|4v)|¬G,|;{A4}|¬G{H11}({H9}(u|4|↔n|
13520  4v)|4α_↓|4(u/v){H11}){H9}/(u/v)|¬G.|;{A9}{H9L11M29}|π|9|1|≡8
13521  |≡.|4|9[|ε|*/M|↔P|↔l|\|π] Estimate the accuracy 
13525  of the ``improved'' double-precision subroutines 
13530  of exercises 4 and 5, in the sense of exercise 
13540  7.|'{A3}|H*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/
folio 328 galley 17
13542  A.W.)!f. 328!ch. 4!g. 19a|'{A12}!|4|4|4|4The 
13547  di=erence between exponents had a behavior approximately 
13554  as the probabilities given in the following table, 
13562  for various radices |εb:|'{A12}{H9L11M29}{M22.6}|∂!!!!!!|∂!!
13566  !!!!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;|>|¬Ge|βu|4α_↓|4e|βv|¬G|;
13569  b|4α=↓|42|;b|4α=↓|410|;b|4α=↓|416|;b|4α=↓|464|;
13573  >{B2}|J#>|>0|;0.33|;0.47|;0.47|;0.56|;>|>1|;0.12|;
13585  0.23|;0.26|;0.27|;>|>2|;0.09|;0.11|;0.10|;0.04|;
13595  >|>3|;0.07|;0.03|;0.02|;0.02|;>|>4|;0.07|;0.01|;
13607  0.01|;0.02|;>|>5|;0.04|;0.01|;0.02|;0.00|;>|>
13618  |πover|45|;0.28|;0.13|;0.11|;0.09|;>{B2}|J#>|>
13626  average|;3.1|9|1|;0.9|9|1|;0.8|9|1|;0.5|9|1|;
13631  >{A12}{H10L12M29}(The ``over 5'' line includes 
13637  essentially all of the cases when one operand 
13645  is zero, but the ``average'' does not include 
13653  these cases.)|'!|4|4|4|4When |εu |πand |εv |πhave 
13660  the same sign and are normalized, then |εu|4α+↓|4v 
13668  |πeither requires one shift to the |εright |π(for 
13676  fraction over⊗ow), or no normalization shifts 
13682  whatever. When |εu |πand |εv |πhave opposite 
13689  signs, we have zero or more |εleft |πduring the 
13698  normalization. The following table gives the 
13704  observed number of shifts required:|'{A12}{H9L11M24}|∂!!!!!!
13709  !!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂!!!!!!|∂|E|;|>|;|εb|4α=↓|42|;
13713  b|4α=↓|410|;b|4α=↓|416|;b|4α=↓|464|;>{B2}|J#>
13718  |>|9|πShift|4right|41|'0.20|;0.07|;0.06|;0.02|;
13724  >|>|9No|4shift|'0.59|;0.80|;0.82|;0.87|;>|>|9Shift|4left|41|
13733  '0.07|;0.08|;0.07|;0.06|;>|>|9Shift|4left|42|'
13741  0.03|;0.02|;0.01|;0.01|;>|>|9Shift|4left|43|'
13748  0.02|'0.00|;0.01|;0.00|;>|>|9Shift|4left|44|'
13755  0.02|'0.01|;0.00|;0.01|;>|>|9Shift|4left|4|→|¬Q4|'
13762  0.06|;0.02|;0.02|;0.02|;>t{A12}{H10L12M29}The 
13768  last line includes all cases where the result 
13776  was zero. The average number of left shifts per 
13785  normalization was about 0.9 when |εb|4α=↓|4|π2; 
13791  0.2 when |εb|4α=↓|410 |πor 16; and 0.1 when |εb|4α=↓|464.|'
13800  {A12}|π|≡B|≡.|9|≡T|≡h|≡e |≡f|≡r|≡a|≡c|≡t|≡i|≡o|≡n 
13802  |≡p|≡a|≡r|≡t|≡s|≡.|9|4Further analysis of ⊗oating-point 
13806  routines can be based on the |εstatistical distribution 
13814  of the fraction parts |πof randomly chosen normalized 
13822  ⊗oating-point numbers. In this case the facts 
13829  are quite surprising, and there is an interesting 
13837  theory which accounts for the unusual phenomina 
13844  which are observed.|'!|4|4|4|4For convenience 
13849  let us temporarily assume that we are dealing 
13857  with ⊗oating-|εdecimal |π(i.e., radix 10) arithmetic; 
13863  modi_cations of the following discussion to any 
13870  other positive integer base |εb |πwill be very 
13878  straightforward. Suppose we are given a ``random'' 
13885  positive normalized number (|εe,|4f)|4α=↓|410|ge|4|¬O|4f. 
13889  |πSince |εf |πis normalized, we know that its 
13897  leading digit is 1, 2, 3, 4, 5, 6, 7, 8, or 9, 
13910  and it seems natural to assume that each of these 
13920  leading digits will occur about one-ninth of 
13927  the time. But, in fact, the behavior in practice 
13936  is quite di=erent. For example, the leading digit 
13944  tends to be equal to 1 over 30 percent of the 
13955  time*3|'!|4|4|4|4One way to test the assertion 
13962  just made is to take a table of physical constants 
13972  (e.g., the speed of light, the force of gravity) 
13981  from some standard reference. If we look at the 
13990  |εHandbook of Mathematical Functions |π(U.S. 
13995  Dept of Commerce, 1964), for example, we _nd 
14003  that 8 of the 28 di=erent physical constants 
14011  given in Table 2.3, roughly 29 percent, have 
14019  leading digit equal to 1. The decimal values 
14027  of |εn*3 |πfor |ε1|4|¬E|4n|4|¬E|4100 |πinclude 
14032  exactly 30 entries beginning with 1; so do the 
14041  decimal values of 2|g|εn |πand of |εF|βn, for 
14049  |ε1|4|¬E|4n|4|¬E|4100. |πWe might also try looking 
14055  at census reports, or a Farmer's Almanack, etc..n(bbbbbb*?but
14062   not a telephone directory).|'!|4|4|4|4In the 
14069  days before pocket calculators, the pages in 
14076  well-used tables of logarithms tended to get 
14083  quite dirty in the front, while the last pages 
14092  stayed relatively clean and neat. This phenomenon 
14099  was apparently _rst mentioned by the American 
14106  astronomer Simon Newcomb [|εAmer. J. Math. |≡4 
14113  (1881), 39<40], |πwho gave good grounds for believing 
14121  that the leading digit |εd |πoccurs with probability 
14129  log|β1|β0(1|4α+↓|41/|εd). |πThe same distribution 
14133  was discovered empirically, many years later, 
14139  by Frank Benford, who reported the results of 
14147  20,229 observations taken from di=erent sources 
14153  [|εProc. Amer. Philosophical Soc. |≡7|≡8 (1938), 
14159  551<572].|'|π!|4|4|4|4In order to account for 
14165  this leading-digit law, let's take a closer look 
14173  at the way we write numbers in ⊗oating-point 
14181  notation. If we take any positive number |εu, 
14189  |πits leading digits are determined by the value 
14197  (log|β1|β0|4|εu) |πmod 1: The leading digit is 
14204  less than |εd |πif and only if|'{A9}(log|β1|β0|4|εu)|πmod|41
14211  |4|¬W|4log|β1|β0|4|εd,|J!(1)|;{A9}|πsince 10|εf|βu|4α=↓|410|
14213  ur(|πlog|β1|β0|4|εu)|πmod|41|)|).|'!|4|4|4|4Now 
14215  if we have any ``random'' positive number |εU, 
14223  |πchosen from sone reasonable distribution such 
14229  as might occur in nature, we might expect that 
14238  (log|β1|β0|4|εU)|πmod|41 would be uniformly distributed 
14243  between zero and one, at least to a very good 
14253  approximation. (Similarly, we expect |εU |πmod 
14259  1, |εU|g2 |πmod 1, |ε|¬K|v4U|4α+↓|4|≤p|) |πmod 
14265  1, etc., to be uniformly distributed. We expect 
14273  a roulette wheel to be unbiased, for essentially 
14281  the same reason.) Therefore by (1) the leading 
14289  digit will be 1 with probability log|β1|β0|42|4|¬V|430.103 
14296  percent; it will be 2 with probability log|β1|β0|43|4α_↓|4lo
14303  g|β1|β0|42|4|¬V|417.609 percent; and, in general, 
14308  if |εr |πis any real value between 1 and 10, 
14318  we ought to have |ε10f|βU|4|¬E|4r |πapproximately 
14324  log|β1|β0|4|εr |πof the time.|'!|4|4|4|4Another 
14329  way to explain this law is to say that a random 
14340  value |εU |πshould appear at a random point on 
14349  a slide rule, according to the uniform distribution, 
14357  since the distance from the left end of a slide 
14367  rule to the position of |εU |πis proportional 
14375  to (log|β1|β0|4|εU)|πmod 1. The analogy between 
14381  slide rules and ⊗oating-point calculation is 
14387  very close when multiplication and division are 
14394  being considered.|'!|4|4|4|4The fact that leading 
14400  digits tend to be small is important to keep 
14409  in mind; it makes the most obvious techniques 
14417  of ``average error'' estimation for ⊗oating-point 
14423  calculations invalid. The relative error due 
14429  to rounding is usually a little more than expected.|'
14438  !|4|4|4|4Of course, it may justly be said that 
14446  the heuristic argument above does not prove the 
14454  stated law. It merely shows us a plausible reason 
14463  why the leading digits behave the way they do. 
14472  An interesting approach to the analysis of leading 
14480  digits has been suggested by R. Hamming: Let 
14488  |εp(r) |πbe the probability that |ε10f|βU|4|¬E|4r, 
14494  |πwhere |ε1|4|¬E|4r|4|¬E|4|¬E|410 |πand |εf|βU 
14498  |πis the normalized fraction part of a random 
14506  normalized ⊗oating-point number |εU. |πIf we 
14512  think of random quantities in the real world, 
14520  we observe that they are measured in terms of 
14529  arbitrary units; and if we were to change the 
14538  de_nition of a meter or a gram, many of the fundamental 
14549  physical constants would have di=erent values. 
14555  Suppose then that all of the numbers in the universe 
14565  are suddenly multiplied by a constant factor 
14572  |εc; |πour universe of random ⊗oating-point quantities 
14579  should be essentially unchanged by this transformation, 
14586  so{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W)!ch.
folio 330 galley 18
14587   4!f. 330!g. 18a|'{A12}!|4|4|4|4Multiplying everything 
14593  by |εc |πchanges (log|β1|β0|4|εU)|πmod 1 to (log|β1|β0|4|εU|
14599  4α+↓|4|πlog|β1|β0|4|εc)|πmod|41. It is now time 
14604  to set up the formulas which describe the desired 
14613  behavior; we may assume that |ε1|4|¬E|4c|4|¬E|410. 
14619  |πBy de_nition,|'{A9}|εp(r)|4α=↓|4|πprobability{H12}({H10}(l
14621  og|β1|β0|4|εU)|πmod|41|4|¬E|4log|β1|β0|4|εr{H12}){H10}.|;
14622  {A9}|πBy our assumption, we should also have|'
14629  {A9}|εp(r)|4|∂α=↓|4|πprobability{H12}({H10}(log|β1|β0|4|εU|4
14629  α+↓|4|πlog|β1|β0|4|εc)|πmod|41|4|¬E|4log|β1|β0|4|εr{H12})|'
14630  {A4}{A20}|Lα=↓|E>{B20}|L!!|πprobability{H12}({H10}(log|β1|β0
14631  |4|εU)|πmod|41|4|¬E|4log|β1|β0|4|εr|4α_↓|4|πlog|β1|β0|4|εc>
14632  |πor!!(log|β1|β0|4|εU)|πmod|41|4|¬R|41|4α_↓|4log|β1|β0|4|εc{
14632  H12}){H10},!!|πif!!|εc|4|¬E|4r;|?{A4}|L!!|πprobability{H12}(
14633  {H10}1|4α_↓|4log|β1|β0|4|εc|4|¬E|4(|πlog|β1|β0|4|εU)|πmod|41
14633  |4|¬E|41|4α+↓|4log|β1|β0|4|εr|4α_↓|4|πlog|β1|β0|4|εc{H12}){H
14633  10},>|πif!!|εc|4|¬R|4r;|?{A6}{A8}|Lα=↓|E>{B8}|L!!p(r/c)|4α+↓
14636  |41|4α_↓|4p(10/c),!!|πif!!|εc|4|¬E|4r;>{A4}|L!!p(10r/c)|4α_↓
14637  |4p(10/c),!!!|4|1|1|1|πif!!|εc|4|¬R|4r.|J!(2)>
14638  {A9}|πLet us now extend the function |εp(r) |πto 
14646  values outside the range |ε1|4|¬E|4r|4|¬E|410, 
14651  |πby de_ning |εp(10|gnr)|4α=↓|4p(r)|4α+↓|4n; 
14654  |πthen if we replace 10/|εc |πby |εd, |πthe last 
14663  equation of (2) may be written|'{A9}|εp(rd)|4α=↓|4p(r)|4α+↓|
14669  4p(d).|J!(3)|;{A9}|π{H10L12M29}If our assumption 
14673  about invariance of the distribution under multiplication 
14680  by a constant factor is valid, then Eq. (3) must 
14690  hold for all |εr|4|¬Q|40 |πand |ε1|4|¬E|4d|4|¬E|410. 
14696  |πThe facts that |εp(1)|4α=↓|40, p(10)|4α=↓|41 
14701  |πnow imply that|'{A9}|ε1|4α=↓|4p(10)|4α=↓|4p{H12}({H10}(|¬K
14704  |v210|))|gn{H12}){H10}|4α=↓|4p(|¬K|v210|))|4α+↓|4p{H12}({H10
14704  }(|¬K|v210|))|gn|gα_↓|g1|4α=↓|1|1|¬O|4|¬O|4|¬O|1|1α=↓|4np(|¬
14704  K|v210|));|;{A9}|πhence we deduce that |εp(10|gm|g/|gn)|4α=↓
14709  |4m/n |πfor all positive integers |εm |πand |εn. 
14717  |πIf we now decide to require that |εp |πis continuous, 
14727  we are forced to conclude that |εp(r)|4α=↓|4log|β1|β0|4r, 
14734  |πand this is the desired law.|'{H10L12M29}!|4|4|4|4Although
14740   this argument may be more convincing than the 
14749  _rst one, it doesn't really hold up under scrutiny 
14758  if we stick to conventional notions of probability. 
14766  The traditional way to make the above argument 
14774  rigorous is to assume that there is some underlying 
14783  distribution of numbers |εF(u) |πsuch that a 
14790  given positive number |εU |πis |→|¬E|εu |πwith 
14797  probability |εF(u); |πand that|'{A9}|εp(r)|4α=↓|4|↔K|uc|)m|)
14801  |4{H12}({H10}F(10|gmr)|4α_↓|4F(10|gm){H12}){H10}|J!(4)|;
14802  {A9}summed over all values |ε|→α_↓|¬X|4|¬W|4m|4|¬W|4|¬X. 
14807  |πOur assujptions about scale invariance and 
14813  continuity led to the conclusion that|'{A9}|εp(r)|4α=↓|4|πlo
14819  g|β1|β0|4|εr.|;{A9}|πUsing the same argument, 
14824  we could ``prove'' that|'{A9}|ε|↔K|uc|)m|)|4{H12}({H10}F(b|g
14828  mr)|4α_↓|4F(b|gm){H12}){H10}|4α=↓|4|πlog|ε|βb|4r,|J!(5)|;
14829  {A9}{H10L12M29}|πfor each integer |εb|4|¬R|42, 
14833  |πwhen |ε1|4|¬E|4r|4|¬E|4b. |πBut there |εis 
14838  |πno distribution function |εF |πwhich satis_es 
14844  this equation for all such |εb |πand |εr*3|'|π!|4|4|4|4One 
14853  way out of the di∃culty is to regard the logarithm 
14863  law |εp(r)|4α=↓|4|πlog|β1|β0|4|εr |πas only a 
14868  very close |εapproximation |πto the true distribution. 
14875  The true distribution itself may perhaps be changing 
14883  as the universe expands, becoming a better and 
14891  better approximation as time goes on; and if 
14899  we replace 10 by an arbitrary base |εb, |πthe 
14908  approximation appears to be less accurate (at 
14915  any given time) as |εb |πgets larger. Another 
14923  rather appealing way to resolve the dilemma, 
14930  by abandoning the traditional idea of a distribution 
14938  function, has been suggested by R. A. Raimi, 
14946  |εAMM |≡7|≡6 (1969), 342<348.|'|π!|4|4|4|4The 
14951  hedging in the last paragraph is probably a very 
14960  unsatisfactory explanation, and so the following 
14966  further calculation (which sticks to regorous 
14972  mathematics and avoids any intuitive, yet paradoxical, 
14979  notions of probability) should be welcome. Let 
14986  us consider the distribution of the leading digits 
14994  of the |εpositive integers, |πinstead of the 
15001  distribution for some imagined set of real numbers. 
15009  The investigation of this topic is quite interesting, 
15017  not only because it sheds some light on the probability 
15027  distributions of ⊗oating-point data, but also 
15033  because it makes a particularly instructive example 
15040  of how to combine the methods of discrete mathematics 
15049  with the methods of in_nitesimal calculus.|'!|4|4|4|4In 
15056  the following discussion, let |εr |πbe a _xed 
15064  real number, |ε1|4|¬E|4r|4|¬E|410; |πwe will 
15069  atempt to make a reasonable attempt at de_ning 
15077  |εp(r), |πthe ``probability'' that the representation 
15083  10|ε|ge|rN|4|¬O|4f|βN |πof a ``random'' positive 
15088  integer |εN |πhas |ε10f|βN|4|¬W|4r, |πassuming 
15093  in_nite precision.|'!|4|4|4|4To start, let us 
15099  try to _nd the probability using a limiting method 
15108  like the de_nition of ``Pr'' in Section 3.5. 
15116  One nice way to rephrase that de_nition is to 
15125  de_ne|'{A9}|h|ε!|9P|β0(n)|4α=↓|4!|∂oikjuyh|E|n|'
15127  |L1,!!|πif|4|1|εn|4α=↓|410|ge|4|¬O|4f|4|1|πwhere|4|110|εf|4|
15127  ¬W|4r;>!|9P|β0(n)|4α=↓|E|'|πi.e.,|4|1if|4|1(log|β1|β0|4|εn)|
15129  πmod|41|4|¬W|4log|β1|β0|4|εr;!!(6)|?|L0,!!|πotherwise.>
15131  {A9}|πNow |εP|β0(1), P|β0(2),|4.|4.|4. |πis an 
15136  in_nite sequence of zeros and ones, with ones 
15144  to represent the cases which contribute to the 
15152  probability. We can try to ``average out'' this 
15160  sequence, by de_ning|'{A9}|εP|β1(n)|4α=↓|4|(1|d2n|)|4|↔K|uc|
15163  )1|¬Ek|¬En|)|4P|β0(k).|J!(7)|;{A9}|πThus if we 
15167  generate a random integer between 1 and |εn |πusing 
15176  the techniques of Chapter 3, and convert it to 
15185  ⊗oating-decimal form (|εe,|4f), |πthe probability 
15190  that |ε10f|4|¬W|4r |πis exactly |εP|β1(n). |πIt 
15196  is natural to let lim|ε|βn|β|¬M|β|¬X P|β1(n) 
15202  |πbe the ``probability'' |εp(r) |πwe are after, 
15209  and that is just what we did in Section 3.5.|'
15219  !|4|4|4|4But in this case the limit does not 
15227  exist: For example, let us consider the subsequence|'
15235  {A9}|εP|β1(s), P|β1(10s), P|β1(100s),|4.|4.|4.|4,|4P|β1(10|g
15237  ns),|4.|4.|4.|4,|;{U0}{H10L12M29}58320#Computer 
folio 333 galley 19
15239  Programming!(Knuth/A.W.)!F. 333!ch. 4!g. 19a|'
15243  {A12}|πwhere |εs |πis  a real number, |ε1|4|¬E|4s|4|¬E|410. 
15251  |πIf |εs|4|¬E|4r, |πwe _nd that|'{A9}|εP|β1(10|gns)|'
15257  {A4}α=↓|4|(1|d210|gns|)|4(|"pr|"P|4α_↓|41|4α+↓|4|"p10r|"P|4α
15257  _↓|410|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4|"p10|gn|gα_↓|g1r|"P|4α_↓|4
15257  10|gn|gα_↓|g1|4α+↓|4|"l10|gns|"L|'{A3}α+↓|41|4α_↓|410|gn)|?
15259  {A4}α=↓|4|(1|d210|gns|)|4{H12}({H10}r(1|4α+↓|410|4α+↓|4|¬O|4
15259  |¬O|4|¬O|4α+↓|410|gn|gα_↓|g1)|4α+↓|4O(n)|4α+↓|4|"l10|gns|"L|
15259  4α_↓|41|4α_↓|410|4α_↓|1|1|¬O|4|¬O|4|¬O|1|1α_↓|410|gn{H12})|'
15260  {A4}{H10}α=↓|4|(1|d210|gns|)|4{H12}({H10}|f1|d39|)(10|gnr|4α
15260  _↓|410|gn|gα+↓|g1)|4α+↓|4|"l10|gns|"L|4α+↓|4O(n){H12}){H10},
15260  |J!(8)|'>{A9}|πwhere |εr|4α=↓|4r|β0|4|¬O|4r|β1r|β2|4.|4.|4. 
15264  |πin decimal notation. As |εn|4|¬M|4|¬X, P|β1(10|gns) 
15270  |πtherefore approaches the limiting value |ε1|4α+↓|4(r|4α_↓|
15275  410)/9s. |πThe above calculation for the case 
15282  |εs|4|¬E|4r |πcan be modi_ed so that it is valid 
15291  for |εs|4|¬Q|4r |πif we replace |ε|"l10|gns|"L|4α+↓|41 
15297  |πby |"p10|ε|gnr|"P; |πso when |εs|4|¬R|4r, |πwe 
15303  _nd that the limiting value is |ε10(r|4α_↓|41)/9s. 
15310  [|πSee J. Franel, |εNaturforschende Gesellschaft, 
15315  Vierteljahrsschrift |≡6|≡2 (|πZ|=4urich, 1917), 
15319  286<295.]|'{H10L12M29}!|4|4|4|4Therefore |εP|β1(n) 
15322  |πhas subsequences whose limit goes from |ε(r|4α_↓|41)/9 
15329  |πup to |ε10(r|4α_↓|41)/9r |πand down again to 
15336  |ε(r|4α_↓|41)/9, |πas |εs |πgoes from 1 to |εr 
15344  |πto 10. We see that |εP|β1(n) |πhas no limit, 
15353  and |εP|β1(n) |πis not a particularly good approximation 
15361  to ur conjectured answer log|β1|β0|4|εr |πeither*3|'
15367  !|4|4|4|4Since |εP|β1(n) |πdoesn't approach a 
15372  limit, we can try to use the same idea as (7) 
15383  once again, to ``average out'' this anomalous 
15390  behavior. In general, let|'{A9}|εP|βm|βα+↓|β1(n)|4α=↓|4|(1|d
15394  2n|)|4|↔K|uc|)1|¬Ek|¬En|)|4P|βm(k).|J!(9)|;{A9}{H10L12M29}|π
15395  Then |εP|βm|βα+↓|β1(n) |πwill tend to be a more 
15403  well-behaved sequence than |εP|βm(n). |πLet us 
15409  try to examine the behavior of |εP|βm|βα+↓|β1(n), 
15416  |πfor large |εn. |πOur experience above with 
15423  the special case |εm|4α=↓|40 |πindicates that 
15429  it might be worth while to consider the subsequence 
15438  |εP|βm|βα+↓|β1(10|gns); |πand the following results 
15443  can now be derived:|'{A12}|≡L|≡e|≡m|≡m|≡a |≡Q|≡.|9|4|εFor 
15449  any integer m|4|¬R|41 and any real number |≤e|4|¬Q|40, 
15457  there are functions Q|βm(s), R|βm(s) and an integer 
15465  N|βm(|≤e), such that whenever n|4|¬Q|4N|βm(|≤e) 
15470  and 1|4|¬E|4s|4|¬E|410, we have|'{A17}(10)|E|?
15475  {B8}|¬GP|βm(10|gns)|4α_↓|4Q|βm(s)|¬G|4|¬W|4|≤e,!!|πif!!|εs|4
15475  |¬E|4r;|;{A4}|¬GP|βm(10|gns)|4α_↓|4{H12}({H10}Q|βm(s)|4α+↓|4
15476  R|βm(s){H12}){H10}|¬G|4|¬W|4|≤e,!!|πif!!|εs|4|¬Q|4r.|;
15477  {A9}Furthermore the functions Q|βm(s), R|βm(s) 
15482  satisfy the relations|'{A9}Q|βm(s)|4|∂α=↓|4|(1|d2s|)|4|↔a|(1
15485  |d29|)|4|↔j|ur10|)1|)Q|βm|βα_↓|β1(t)|4dt|4α+↓|4|↔j|urs|)1|)Q
15485  |βm|βα_↓|β1(t)|4dt|4α+↓|4|(1|d29|)|4|↔j|ur10|)r|)|4R|βm|βα_↓
15485  |β1(t)|4dt|↔s;|;{A4}| R|βm(s)|4|Lα=↓|4|(1|d2s|)|4|↔j|urs|)r|
15486  )|4R|βm|βα_↓|β1(t)|4dt;|J!(11)>{A4}| Q|β0(s)|4|Lα=↓|41,!!R|β
15487  0(s)|4α=↓|4|→α_↓1.>{A12}Proof.|9|4|πConsider 
15489  the functiohs |εQ|βm(s), R|βm(s) |πde_ned by 
15495  (11), and let|'{A9}|h|ε|∂S|βm(t)|4α=↓!|∂Q|βm(t)|4α+↓|4R|βm(t
15498  ),!!|∂t|4|¬E|4r.|E|n|;{A8}|LS|βm(t)|4α=↓|E>{B8}|L|LQ|βm(t),|
15500  Lt|4|¬E|4r.>{A4}|L|LQ|βm(t)|4α+↓|4R|βm(t),!!t|4|¬Q|4r.>
15502  {A9}{H10L12M29}|πWe will prove the lemma by induction 
15509  on |εm.|'{H10L12M29}|π!|4|4|4|4First, let |εm|4α=↓|41; 
15514  |πthen |εQ|β1(s)|4α=↓|4{H12}({H10}1|4α+↓|4(s|4α_↓|41)|4α+↓|4
15515  (r|4α_↓|410)/9{H12}){H10}/s|4α=↓|41|4α+↓|4(r|4α_↓|410)/9s, 
15516  |πand |εR|β1(s)|4α=↓|4(r|4α_↓|4s)/s. |πFrom (8) 
15520  we _nd that|'{H10}|ε|¬GP|β1(10|gns)|4α_↓|4S|β1(s)|¬G|4α=↓|4O
15523  (n)/10|gn;|;{A9}|πthis establishes the lemma 
15528  when |εm|4α=↓|41.|'|π!|4|4|4|4Now for |εm|4|¬Q|41, 
15533  |πwe have |εP|βm(10|gns)|4α=↓|'{A6}|(1|d2s|)|4|↔a|↔k|uc|)0|¬
15536  Ej|¬Wn|)|4|(1|d210|gn|gα_↓|gj|)|4|↔k|uc|)10|gj|¬Ek|¬W10|gj|g
15536  α+↓|g1|)|4|(1|d210|gj|)|4P|βm|βα_↓|β1(k)|4α+↓|4|↔k|uc|)10|gn
15536  |¬Ek|¬E10|gns|)|4|(1|d210|gn|)|4P|βm|βα_↓|β1(k)|↔s,|;
15537  {A9}|πand we want to approximate this quantity. 
15544  By induction the di=erence|'{A9}|ε|↔g|4|↔k|uc|)10|gj|¬Ek|¬E1
15548  0|gjq|)|4|(1|d210|gj|)|4P|βm|βα_↓|β1(k)|4α_↓|4|↔k|uc|)10|gj|
15548  ¬Ek|¬E10|gjq|)|4|(1|d210|gj|)|4S|βm|βα_↓|β1|↔a|(k|d210|gj|)|
15548  ↔s|↔g|J!(13)|;{A9}|πis less than |εq|≤e |πwhen 
15554  |ε1|4|¬E|4q|4|¬E|410 |πand |εj|4|¬Q|4N|βm|βα_↓|β1(|≤e); 
15557  |πand since |εS|βm|βα_↓|β1(t) |πis continuous, 
15562  it is a Riemann-integrable function, and the 
15569  di=erence|'{A9}|ε|↔g|4|↔k|uc|)10|gj|¬Ek|¬E10|gjq|)|4|(1|d210
15570  |gj|)|4S|βm|βα_↓|β1|↔a|(k|d210|gj|)|↔s|4α_↓|4|↔j|urq|)1|)|4S
15570  |βm|βα_↓|β1(t)|4dt|4|↔g|J!(14)|;{A9}|πis less 
15573  than |ε|≤e |πfor all |εj |πgreater than some 
15581  number |εN, |πindependent of |εq, |πby the de_nition 
15589  of integration. We may choose |εN |πto be |→|¬Q|εN|βm|βα_↓|β
15597  1(|≤e). |πTherefore for |εn|4|¬Q|4N, |πthe di=erence|'
15603  {A9}|ε|↔g|4P|βm(10|gns)|4α_↓|4|(1|d2s|)|4|↔a|↔k|uc|)0|¬Ej|¬W
15603  n|)|4|(1|d210|gn|gα_↓|gj|)|4|↔j|ur10|)1|)|4S|βm|βα_↓|β1(t)|4
15603  dt|4α+↓|4|↔j|urs|)1|)|4S|βm|βα_↓|β1(t)|4dt|↔s|↔g|J!(15)|;
15604  {A9}|πis bounded by|'{A9}|ε|↔k|uc|)0|¬Ej|¬EN|)|4|(M|d210|gn|
15607  gα_↓|gj|)|4α+↓|4|↔k|uc|)N|¬Wj|¬Wn|)|4|(11|≤e|d210|gn|gα_↓|gj
15607  |)|4α+↓|411|≤e,|;{A9}|πif |εM |πis an upper bound 
15614  for (13)|4α+↓|4(14) which is valid for all |εj. 
15622  |πFinally, the sum |¬K|ur|)|ε0|¬Ej|¬Wn|)|4(1/10|gn|gα_↓|gj),
15625   |πwhich appears in (15), is equal to |ε(1|4α_↓|41/10|gn)/9;
15633   |πso|'{A9}|ε|↔g|4P|βm(10|gns)|4α_↓|4|(1|d2s|)|4|↔a|(1|d29|)
15635  |4|↔j|ur10|)1|)S|βm|βα_↓|β1(t)|4dt|4α+↓|4|↔j|urs|)1|)S|βm|βα
15635  _↓|β1(t)|4dt|↔s|↔g|;{A9}{H10}|πcan be made smaller 
15640  than, say, 20|ε|≤e, |πif |εn |πis taken large 
15648  enough. comparing this with (10) and (11) completes 
15656  the proof.|'{A12}!|4|4|4|4The gist of Lemma Q 
15663  is that we have the limiting relationship|'{A9}|uwlim|uc|)|ε
15670  n|¬M|¬X|)|4P|βm(10|gns)|4α=↓|4S|βm(s).|J!(16)|;
15671  {A9}|H*?{U0}{H10L12M29}58320#Computer Programming!(Knuth/A.-W
folio 336 galley 20
15672  .)!f. 336!ch. 4!g. 20a|'{A12}|πAlso, since |εS|βm(s) 
15679  |πis not constant as |εs |πvaries, the limit|'
15687  {A9}|uwlim|uc|)|εn|¬M|¬X|)|4P|βm(n),|;{A9}|πwhich 
15689  would be our desired ``probability,'' does not 
15696  exist for any |εm. |πThe situation is shown in 
15705  Fig. 5, which shows the values of |εS|βm(s) |πwhen 
15714  |εm |πis small and |εr|4α=↓|42.|'|π!|4|4|4|4Even 
15720  though |εS|βm(s) |πis not a constant, so that 
15728  we do not have a de_nite limit for |εP|βm(n), 
15737  |πnote that already for |εm|4α=↓|43 |πin Fig. 
15744  5 the value of |εS|βm(s) |πstays very close to 
15753  log|β1|β0|42|4α=↓|40.30103|4.|4.|4.|4. Therefore 
15755  we have good reason to suspect that |εS|βm(s) 
15763  |πis very close to log|β1|β0|4|εr |πfor all large 
15771  |εm, |πand, in fact, that the sequence of functions 
15780  |ε|↔bS|βm(s)|↔v |πconverges uniformly to the 
15785  constant function log|β1|β0|4|εr.|'|π!|4|4|4|4It 
15789  is interesting to prove this conjecture by explicitly 
15797  calculating |εQ|βm(s) Q|βm(s) |πand |εR|βm(s) 
15802  |πfor all |εm, |πas in the proof of the following 
15812  theorem:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡F|≡.|9|4|εLet 
15815  S|βm(s) be the limit de⊂ned in (16). For all 
15824  |≤e|4|¬Q|40, there exists a number N(|≤e) such 
15831  that|'{A9}|¬GS|βm(s)|4α_↓|4|πlog|β1|β0|4|εr|¬G|4|¬W|4|≤e!!|π
15832  for!!|ε1|4|¬E|4s|4|¬E|410,|J!(17)|;{A9}whenever 
15834  m|4|¬Q|4N(|≤e).|'{A6}{H9L11M29}|π|≡F|≡i|≡g|≡. 
15836  |≡5|≡.|9|4The probability that the leading digit 
15842  is 1.|;{A6}{H10L12M29}|εProof.|4|9|πIn view of 
15847  Lemma Q, we can prove this result if we can show 
15858  there is a number |εM |πdepending on |ε|≤e |πsuch 
15867  that, for |ε1|4|¬E|4s|4|¬E|410,|'{A9}|¬GQ|βm(s)|4α_↓|4|πlog|
15870  β1|β0|4|εr|¬G|4|¬W|4|≤e!!|πand!!|ε|¬GR|βm(s)|¬G|4|¬W|4|≤e|J!
15870  (18)|;{A9}|πfor all |εm|4|¬Q|4M.|'!|4|4|4|4|πIt 
15875  is not di∃cult to solve the recurrence formula 
15883  (11) for |εR|βm: |πWe have |εR|β0(s)|4α=↓|4|→α_↓1, 
15889  R|β1(s)|4α=↓|4|→α_↓1|4α+↓|4r/s, R|β2(s)|4α=↓|4|→α_↓1|4α+↓|4(
15890  r/s){H12}({H10}1|4α+↓|4|πln|4(|εs/r){H12}){H10}, 
15891  |πand in general|'{A9}|εR|βm(s)|4α=↓|4|→α_↓1|4α+↓|4|(r|d2s|)
15894  |4{H12}|↔a{H10}1|4α+↓|4|(1|d21*3|)|4|πln|4|↔a|(|εs|d2r|)|↔s|4
15894  α+↓|4|(1|d2(m|4α_↓|41)*3|)|4|↔a|πln|4|↔a|(s|d2r|)|↔s|↔s|gm|gα
15894  _↓|g1{H12}|↔s{H10}.|J!(19)|;{A9}|πFor the stated 
15898  range of |εs, |πthis converges uniformly to|'
15905  {A9}|ε|→α_↓1|4α+↓|4(r/s)|4|πexp|4{H12}({H10}ln|4(|εs/r){H12}
15905  )|4{H10}α=↓|40.|;{A9}{H10L12M29}|π!|4|4|4|4The 
15907  recurrence (11) for |εQ|βm |πtakes the form|'
15914  {A9}|εQ|βm(s)|4α=↓|4|(1|d2s|)|4|↔ac|βm|4α+↓|41|4α+↓|4|↔j|urs
15914  |)1|)Q|βm|βα_↓|β1(t)|4dt|↔s,|J!(20)|;{A9}|πwhere|'
15916  {A9}|εc|βm|4α=↓|4|(1|d29|)|4|↔a|↔j|ur10|)1|)Q|βm|βα_↓|β1(t)|
15916  4dt|4α+↓|4|↔j|ur10|)r|)R|βm|βα_↓|β1(t)|4dt|↔s|4α_↓|41.|J!(21
15916  )|;{A9}|πThe solution to recurrence (20) is easily 
15924  found by trying out the _rst few cases and guessing 
15934  at a formula which can be proved by induction; 
15943  we _nd that|'{A9}|εQ|βm(s)|4α=↓|41|4α+↓|4|(1|d2s|)|4|↔ac|βm|
15946  4α+↓|4|(1|d21*3|)|4c|βm|βα_↓|β1|4|πln|4|εs|4α+↓|1|1|¬O|4|¬O|4
15946  |¬O|1|1α+↓|4|(1|d2(m|4α_↓|41)*3|)|4c|β1(|πln|4|εs)|gm|gα_↓|g1
15946  |↔s.|J!(22)|;{A9}|π!|4|4|4|4It remains for us 
15951  to calculate the coe∃cients |εc|βm, |πwhich by 
15958  (19), (21), and (22) satisfy the relations|'{A9}|εc|β1|4α=↓|
15965  4(r|4α_↓|410)/9|;{A9}|ε!c|βm|βα+↓|β1|4α=↓|4|(1|d29|)|4{H12}|
15966  ↔a{H10}c|βm|4|πln|410|4α+↓|4|(1|d22*3|)|4|εc|βm|βα_↓|β1(|πln|
15966  410)|g2|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(|ε1|d2m*3|)|4c|β1(|πl
15966  n|410)|ε|gm|'{A4}α+↓|4r|↔a1|4α+↓|4|(1|d21*3|)|4|πln|4|(|ε10|d
15967  2r|)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d2m*3|)|4|↔a|πln|4|ε|(
15967  10|d2r|)|↔s|gm|↔s|4α_↓|410{H12}|↔s{H10}.|J!(23)|?
15968  {A6}|πThis sequence appears to be very complicated, 
15975  but actually we can analyze it without di∃culty 
15983  with the help of generating functions. Let|'{A9}|εC(z)|4α=↓|
15990  4c|β1z|4α+↓|4c|β2z|g2|4α+↓|4c|β3z|g3|4α+↓|1|1|¬O|4|¬O|4|¬O|4
15990  ;|;{A9}|πthen since 10|ε|gz|4α=↓|41|4α+↓|4z|4|πln|410|4α+↓|4
15993  (1/2*3)|4|εz|4|πln|410)|g2|4α+↓|1|1|¬O|4|¬O|4|¬O|4, 
15994  we deduce that|'{A9}|h|εc|βm|βα+↓|β1|4|∂α=↓|4|∂10|9|1c|βm|βα
15997  +↓|β1|4α+↓|4c|βm|4|πln|410|4α+↓|1|1.|4.|4.|1|1α+↓|4|εm*3|4c|β
15997  1(|πln|410)|gm|9|1|E|n|;|ε| c|βm|βα+↓|β1|4|Lα=↓|4|L|(1|d210|
15998  )|4c|βm|βα+↓|β1|4α+↓|4|(9|d210|)|4c|βm|βα+↓|β1>
15999  {A4}|Lα=↓|4|L|(1|d210|)|4|↔ac|βm|βα+↓|β1|4α+↓|4c|βm|4|πln|41
15999  0|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d2|εm*3|)|4c|β1(|πln|410)
15999  |ε|gm|↔s>{A6}|L|Lα+↓|4|(r|d210|)|4{H12}|↔a{H10}1|4α+↓|1|1|¬O
16000  |4|¬O|4|¬O|1|1α+↓|4|(1|d2m*3|)|4|↔a|πln|4|ε|(10|d2r|)|↔s|gm{H
16000  12}|↔s{H10}|4α_↓|41>{A9}  |πis the coe∃cient 
16006  of |εz|gm|gα+↓|g1 |πin the function|'{A9}|ε|f1|d310|)C(z)10|
16011  gz|4α+↓|4|(r|d210|)|4|↔a|(10|d2r|)|↔s|gz|↔a|(z|d21|4α_↓|4z|)
16011  |↔s|4α_↓|4|(z|d21|4α_↓|4z|)|4.|J!(24)|;{A9}|π{H10L12M29}This
16012   condition holds for all values of |εm, |πso 
16021  (24) must equal |εC(z), |πand we obtain the explicit 
16030  formula|'{A9}|εC(z)|4α=↓|4|(|→α_↓z|d21|4α_↓|4z|)|4|↔a|((10/r
16031  )|gz|gα_↓|g1|4α_↓|41|d210|gz|gα_↓|g1|4α_↓|41|)|↔s.|J!(25)|;
16032  {A9}|π!|4|4|4|4We want to study asymptotic properties 
16038  of the coe∃cients of |εC(z), |πto complete our 
16046  analysis. The large parenthesized quantity in 
16052  (25) approaches ln|4(10/|εr)/|πln|410|4α=↓|41|4α_↓|4log|β1|β
16054  0|4|εr |πas |εz|4|¬M|41, |πso we see that|'{A9}|εC(z)|4α+↓|4
16061  |(1|4α_↓|4|πlog|β1|β0|4|εr|d21|4α_↓|4z|)|4α=↓|4R(z)|J!(26)|;
16062  {A9}|πis an analytic function of the complex 
16069  variable |εz |πin the circle|'{A6}|ε|¬Gz|¬G|4|¬W|4|↔g1|4α+↓|
16074  4|(2|≤pi|d2|πln|410|)|1|1|↔g|1|1.|;{A9}In particular, 
16077  |εR(z) |πconverges for |εz|4α=↓|41, |πso its 
16083  coe∃cients approach zero. This proves that the 
16090  coe∃cients of |εC(z) |πbehave like those of (log|β1|β0|4|εr|
16097  4α_↓|41)/(1|4α_↓|4z), |πthat is,|'{A9}|uwlim|uc|)|εm|¬M|¬X|)
16100  |4c|βm|4α=↓|4|πlog|β1|β0|4|εr|4α_↓|41.|;{A9}|π!|4|4|4|4Final
16101  ly, we may combine this with (22), to show that 
16111  |εQ|βm(s) |πapproaches|'{A9}1|4α+↓|4|(log|β1|β0|4|εr|4α_↓|41
16113  |d2s|)|4|↔a1|4α+↓|4|πln|4|εs|4α+↓|4|(1|d22*3|)|4(|πln|4|εs)|g
16113  2|4α+↓|1|1|¬O|4|¬O|4|¬O|↔s|4α=↓|4|πlog|β1|β0|4|εr|;
16114  {A9}|πuniformly for |ε1|4|¬E|4s|4|¬E|410.|'|Hβ*?*?*?{U0}{H10L12
folio 338 galley 21
16117  M29}58320#Computer Programming!(Knuth/A.-W.)!f. 
16119  338!ch. 4!g. 21a|'{A12}|π!|4|4|4|4The above proofs 
16125  of Lemma Q and Theorem F are slight simpli_cations 
16134  and ampli_cations of methods due to B. J. Flehinger, 
16143  |εAMM |≡7|≡3 (1966), 1056<1061. |πMany authors 
16149  have written about the distribution of initial 
16156  digits, showing that the logarithm lw is a good 
16165  approximation for many undelying distributions; 
16170  see the survey by Ralph A. Raimi, |εAMM |≡8|≡3 
16179  (1976), |πto appear, for a comprehensive review 
16186  of the literature. An appealing new explanation 
16193  of the phenomenon as applied to integers has 
16201  been proposed by D. I. A. Cohen, |εJ. Comb. Theory 
16211  |π(A) |≡2|≡0 (1976), to appear. Another interesting 
16218  (and di=erent) treatment of ⊗oating-point distribution 
16224  has been given by Alan G. Konheim, |εMath. Comp. 
16233  |≡1|≡9 (1965), 143<144. |πFor an approach to 
16240  the de_nition of probability under which the 
16247  logarithm law holds exactly over the integers, 
16254  see exercise 17.|'!|4|4|4|4Therefore we have 
16260  established the logarithmic law for integers 
16266  by direct calculation, at the same time seeing 
16274  that it is an extremely good approximation to 
16282  the average behavior although it is never precisely 
16290  achieved.|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'
16292  {A12}{H9L11M29}|1|9|≡1|≡.|4|9[|ε|*/|↔O|↔L|\|π] 
16293  Given that |εu |πand |εv |πare nonzero ⊗oating-decimal 
16301  numbers |εwith the same sign, |πwhat is the approximate 
16310  probability that fraction over⊗ow occurs during 
16316  the calculation of |εu|4|↔V|4v, |πaccording to 
16322  Sweeney's tables?|'{A3}|9|1|≡2|≡.|9|4[|ε|*/|↔M|↔P|\|π] 
16325  Make further tests of ⊗oating-point addition 
16331  and subtraction, to con_rm or improve on the 
16339  accuracy of Sweeney's tables.|'{A3}|9|1|≡3|≡.|9|4[|ε|*/|↔O|↔C
16343  |\|π] What is the probability that the two leading 
16352  digits of a ⊗oating-decimal number are ``23'', 
16359  according to the logarithm law?|'{A3}{A3}|9|1|≡4|≡.|4|9[|ε|*/
16364  |↔O|↔l|\|π] The text points out that the front 
16372  pages of a well-used table of logarithms get 
16380  dirtier than the back pages do. What if we had 
16390  a |εantilogarithm |πtable instead, i.e., a table 
16397  giving the value of |εx |πwhen log|β1|β0 |εx 
16405  |πis given. Which pages of such a table would 
16414  be the dirtiest?|'{A3}|9|1|≡5|≡.|9|4[|ε|*/M|↔P|↔c|\|π] 
16418  Suppose that |εU |πis a real number uniformly 
16426  distributed in the interval |ε0|4|¬W|4U|4|¬W|41. 
16431  |πWhat is the distribution of the leading digits 
16439  of |εU?|'{A3}|π|1|9|≡6|≡.|4|9[|ε|*/|↔P|↔L|\|π] 
16442  If we have binary computer words containing |εn|4α+↓|41 
16450  |πbits, we might use |εp |πbits for the fraction 
16459  part of ⊗oating-binary numbers, one bit for the 
16467  sign, and |εn|4α_↓|4p |πbits for the exponent. 
16474  This means that the range of values representable, 
16482  ile., the ratio of the largest to the smallest 
16491  positive normalized value, is essentially |ε2|g2|gn|gα_↓|gp.
16496   |πThe same computer word could be used to represent 
16506  ⊗oating |εhexadecimal |πnumbers, i.e., ⊗oating-point 
16511  numbers with radix 16, with |εp|4α+↓|42 |πbits 
16518  for the fraction part {H11}({H9}(|εp|4α+↓|42)/4 
16523  |πhexadecimal digitr{H11}){H9} and |εn|4α_↓|4p|4α_↓|42 
16527  |πbits for the exponent; then the range of values 
16536  would be |ε16|g2|in|iα_↓|ip|i2|4α=↓|42|g2|gn|gα_↓|gp, 
16539  |πthe same as before, and with more bits in the 
16549  fraction part. This may sound as if we are getting 
16559  something for noting, but the normalization condition 
16566  for base 16 is weaker in that there may be up 
16577  to three leading zero bits in the fraction part; 
16586  thus not all of the |εp|4α+↓|42 |πbits are ``suc*?*?*?gni_cant 
16595  bits.''|'!!|1|1On the basis of the logarithmic 
16602  law, what are the probabilities that the fraction 
16610  part of a positive normalized radix 16 ⊗oating-point 
16618  number has exactly 0, 1, 2, and 3 leading zero 
16628  bits? Discuss the desirability of hexadecimal 
16634  versus binary.|'{A3}{A3}|9|1|≡7|≡.|9|4[|ε|*/HM|↔P|↔l|\|π] 
16637  Prove that there is no distribution function 
16644  |εF(u) |πwhich satis_es (5) for each integer 
16651  |εb|4|¬R|42, |πand for all real values |εr |πin 
16659  the range |ε1|4|¬E|4r|4|¬E|4b.|'{A3}|π|9|1|≡8|≡.|4|9[|ε|*/M|↔
16662  P|↔L|\|π] Does (10) hold when |εm|4α=↓|40 |πfor 
16669  suitable |εN|β0(|≤e)?|'{A3}{H9L11M29}|9|1|≡9|≡.|9|4[|ε|*/M|↔P
16671  |↔M|\|π] (P. Diaconis.) Prove that lim|ur|)|εm|¬M|¬X|)|4P|βm
16676  (n)|4α=↓|4P|β0(1) |πfor all _xed |εn.|'{A3}|π|≡1|≡0|≡.|4|9[|
16681  ε|*/HM|↔P|↔l|\|π] The text shows that |εc|βm|4α=↓|4|πlog|β1|β
16686  0 |εr|4α_↓|41|4α+↓|4|≤e|βm, |πwhere |ε|≤e|βm 
16690  |πapproaches zero as |εm|1|1|¬M|1|1|¬X. |πobtain 
16695  the next term in the asymptotic expansion of 
16703  |εc|βm.|'|π{A3}|≡1|≡1|≡.|4|9[|ε|*/M|↔O|↔C|\|π] 
16705  Show that if |εU |πis a random variable which 
16714  is distributed according to the logarithmic law, 
16721  then so is 1/|εU.|'{A3}|π|≡1|≡2|≡.|4|9[|ε|*/M|↔P|↔C|\|π] 
16726  (R. W. Hamming.) The purpose of this exercise 
16734  is to show that the result of ⊗oating-point multiplication 
16743  tends to obey the logarithmic law more perfectly 
16751  than the operands do. Let |εU |πand |εV |πbe 
16760  random, normalized, positive ⊗oating-point numbers, 
16765  whose fraction parts are independently distributed 
16771  with the respective density functions |εf|1(x) 
16777  |πand |εg(x). |πThus, |εf|βu|4|¬E|4r |πand |εf|βv|4|¬E|4s 
16783  |πwith probability |ε|¬J|urr|)1/b|)|1|¬J|urs|)1/b|)|1f|1(x)g
16785  (y)|4dy|4dx, |πfor |ε1/b|4|¬E|4r, s|4|¬E|41. 
16789  |πLet |εh(x) |πbe the density function of the 
16797  fraction part of |εU|4α⊗↓|4V |π(unrounded). De_ne 
16803  the |εabnormality A(f) |πof a density function 
16810  |εf |πas the maximum relative error,|'{A9}|εA(f)|4α=↓|4|uw|π
16816  max|uc|)|ε1/b|¬Ex|¬E1|)|4|↔g|4|(f|1(x)|4α_↓|4l(x)|d2l(x)|)|4
16816  |↔g,|;{A9}|πwhere |εl(x)|4α=↓|41/x|4|πln|4|εb 
16819  |πis the density of the logarithmic distribution. 
16826  Prove that |εA(h)|4|¬E|4|πmin{H11}({H9}A(f), 
16829  A(g){H11}){H9}. [|πIn particular, if either factor 
16835  has logarithmic distribution the product does 
16841  also.]|'{A3}|≡1|≡3|≡.|4|9[|ε|*/M|↔P|↔c|\|π] The 
16844  ⊗oating-point multiplication routine, Algotithm 
16848  4.2.1M, requires zero or one left shifts during 
16856  normalization, depending on whether |εf|βuf|βv|4|¬R|41/b 
16861  |πor not. Assuming that the input operands are 
16869  independently distributed according to the logarithmic 
16875  law, what is the probability that no left shift 
16884  is needed for normalization of the result?|'{A3}|≡1|≡4|≡.|9|
16891  4[|ε|*/HM|↔L|↔c|\|π] Let |εU |πand |εV |πbe random, 
16898  normalized, positive ⊗oating-point numbers whose 
16903  fraction parts are independently distributed 
16908  according to the logarithmic law, and let |εp|βk 
16916  |πbe the probability that the di=erence in their 
16924  exponents is |εk. |πAssuming that the distribution 
16931  of the exponents is independent of the fraction 
16939  parts, give an equation for the probability that 
16947  ``fraction over⊗ow'' occurs during the ⊗oating-point 
16953  addition of |εU|4|↔V|4V, |πin terms of the base 
16961  |εb |πand the quantities |εp|β0, p|β1, p|β2,|4.|4.|4.|4. 
16968  |πCompare this result with exercise 1. (Ignore 
16975  rounding.)|'{A3}|≡1|≡5|≡.|4|9[|ε|*/HM|↔P|↔l|\|π] 
16977  Let |εU, V, p|β0, p|β1,|4.|4.|4. |πbe as in exercise 
16986  14, and assume that radix 10 arithmetic is being 
16995  used. Show that regardless of the values of |εp|β0, 
17004  p|β1, p|β2,|4.|4.|4.|4, |πthe sum |εU|4|↔V|4V 
17009  |πwill |εnot |πobey the logarithmic law exactly, 
17016  and in fact the probability that |εU|4|↔V|4V 
17023  |πhas leading digit 1 is always strictly |εless 
17031  |πthan log|β1|β0|42.|'{A3}|≡1|≡6|≡.|4|9[|ε|*/HM|↔P|↔l|\|π] 
17034  (P. Diaconis.) Let |εP|β0(n) |πbe 0 or 1 for 
17043  each |εn, |πand de_ne ``probabilities'' |εP|βm|βα+↓|β1(n) 
17049  |πby repeated averaging, as in (9). Show that 
17057  if lim|ur|)|εn|¬M|¬X|)|4P|β1(n) |πdoes not exist, 
17062  neither does lim|ur|)|εn|¬M|¬X|)|4P|βm(n) |πfor 
17066  any |εm. [Hint|*/:|\ |πProve that |ε(a|β1|4α+↓|1|1|¬O|4|¬O|4|
17071  ¬O|1|1α+↓|4a|βn)|1|1|¬M|1|10 |πand |εa|βn|βα+↓|β1|4|¬E|4a|βn
17073  |4α+↓|4M/n |πfor some _xed constant. |εM|4|¬Q|40 
17079  |πimplies that |εa|βn|4|¬M|40.]|'{A3}|π|≡1|≡7|≡.|9|4[|ε|*/M|↔
17082  P|↔C|\|π] (R. L. Duncan.) Another way to de_ne 
17090  Pr{H11}({H9}|εS(n){H11}){H9} |πis to evaluate 
17094  lim|ur|)|εn|¬M|¬X|) {H11}({H9}(|¬K|ur|)S(k)|4|πand|4|ε1|¬Ek|
17095  ¬En|)|41/k)/H|βn{H11}){H9}, |πsince it can be 
17100  shown that this ``harmonic probability'' exists 
17106  and equals Pr{H11}({H9}|εS(n){H11}){H9} |πwhenever 
17110  the latter exists according to the de_nition 
17117  in Section 3.5.|'|H*?*?*?*?*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!*!
17120